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1.   INTRODUCTION 

1.1   Background 

A  signal,  represented  mathematically  as  a  function  of  one  or  more 
independent  variables,  conveys  information  about  the  state  or  behavior 
of  a  physical  system.  For  example,  speech  is  represented  as  a  function 
of  time.  Signals  are  said  to  be  continuous-time  (CT)  signals  if  they 
are  defined  over  a  continuum  of  times,  and  discrete-time  (DT)  signals  if 
they  are  defined  only  at  a  discrete  set  of  values  of  the  independent 
variable.  In  addition,  if  the  amplitude  of  a  CT  signal  is  also  con- 
tinuous, then  such  a  signal  is  called  an  analog  signal.  Similarly,  if 
the  amplitude  of  a  DT  signal  is  also  discrete,  such  a  signal  is  called  a 
digital  signal.  Therefore,  digital  signals  are  represented  as  sequences 
of  numbers. 

To  extract  any  meaningful  information  from  a  signal,  it  has  to  be 
processed.  To  that  end,  development  of  techniques  to  process  the  sig- 
nals assumes  great  importance.  Usually,  these  techniques  involve 
transforming  a  signal  to  some  other  advantageous  signal  that  facilitates 
its  analysis.  These  transformation  techniques  are  useful  when  we  want 
to  separate  two  or  more  signals,  or  when  a  particular  component  of  the 
signal  has  to  be  enhanced  in  relation  to  the  others,  or  when  a  desired 
parameter  of  the  signal  has  to  be  estimated.  When  a  transformation 
technique  operates  on  an  analog  signal  and  produces  another  analog 
signal  as  its  output,  it  is  called  analog  signal  processing.   If  the 
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transform  operates  on  a  digital  signal  to  produce  another  digital  signal 
as  its  output,  it  is  called  digital  signal  processing  (DSP). 

DSP  is  employed  in  a  variety  of  fields  of  science  and  technology 
such  as  biomedical  engineering,  acoustics,  sonar,  radar,  seismology, 
speech  communication,  data  communication,  nuclear  science,  spectroscopy 
and  many  others.  As  an  example,  when  a  signal  is  transmitted  over  a 
communication  channel,  it  is  corrupted  in  a  number  of  ways,  by,  e.g., 
channel  distortion,  fading,  and  the  insertion  of  background  noise.  One 
of  the  major  functions  at  the  receiving  end  is  to  compensate  for  all  the 
affecting  disturbances.  DSP  techniques  play  an  important  role  in  such 
an  application. 

DSP  systems  are  especially  attractive  because  they  can  be  realized 
with  flexibility  using  digital  computer,  or  they  can  be  realized  in 
hardware  using  digital  components.  They  can  be  used  to  simulate  analog 
systems,  and  most  importantly,  can  sometimes  be  used  to  realize  other- 
wise impossible  analog  signal  transformations. 

The  Fourier  transform  (FT)  provides  a  very  useful  technique  for  use 
in  CT  signal  processing.  Similarly,  the  discrete  Fourier  transform 
(DFT)  plays  a  vital  role  in  DSP.  Spectrum  analysis,  an  important  com- 
ponent of  DSP,  is  one  area  where  the  DFT  finds  extensive  application. 
But,  for  quite  a  long  time,  no  efficient  means  of  implementing  it, 
either  in  hardware  or  in  software,  was  known.  However,  the  disclosure 
in  1965  by  Cooley  and  Tukey  [1]  of  an  efficient  algorithm  to  compute  the 
DFT  revolutionized  signal  processing.  The  class  of  algorithms  proposed 
by  them  has  come  to  be  known  as  the  fast  Fourier  transform  or  the  FFT. 
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The  FFT  algorithm  reduced  the  computation  time  of  the  DFT  phenomenally, 
facilitating  the  commercial  availability  of  special-purpose  signal 
processing   chips  which  can  operate   at   extremely  high  data   rates. 

Despite    its    tremendous    application,    the   DFT  has   an  unattractive 
feature    in   that    it    transforms   a   real-valued   sequence   also    into   a 
complex-valued   sequence. 

In  1983,  R.  N.  Bracewell  [2]  proposed  an  inherently  real-valued 
transform  called  the  Hartley  transform  (HT).  The  new  transform  has  the 
advantage  that  a  real-valued  signal  always  generates  a  real-valued 
transform  signal.  Secondly,  unlike  the  FT,  the  HT  is  symmetric  that 
is,  both  the  forward  and  inverse  transforms  are  identical.  Thirdly,  the 
HT  is  so  closely  related  to  the  FT,  that  one  can  move  from  the  Fourier 
domain  to  the  Hartley  domain  and  vice  versa  in  a  straightforward  manner. 
Bracewell  also  introduced  the  discrete  Hartley  transform  (DHT),  which 
has  all  the  above  mentioned  advantages  over  the  DFT.  But,  the  most 
important  practical  advantage  of  the  DHT  over  the  DFT  is  that  it  is 
computationally   faster   than  the   DFT. 

The  HT  has  an  interesting  history.  While  working  on  problems 
relating  to  transmission  lines,  in  1942,  Hartley  first  proposed  it  under 
the  name  symmetrical  Fourier  identity  [3],  However,  it  remained  rela- 
tively obscure  for  quite  a  long  time  until  Bracewell  revived  it  in  1983. 
Interestingly,  Zhong  De  Wang,  working  independently,  proposed  an  identi- 
cal transform  in  1981  and  developed  many  of  its  mathematical  properties. 
He   called   it   the  W-Transform    [4-6], 


1.2  Implementation  of    the    Present   Work 

All  the  properties  of  the  ITT  and  the  DHT  are  developed,  analogous 
to  the  properties  of  the  FT  and  the  DFT,  respectively.  Decomposition 
formulas  for  all  the  fast  Hartley  transform  (FHT)  algorithms  are  derived 
and  the  computational  cost  of  different  algorithms  are  compared. 
Studies  in  the  enhancement  of  signal-to-noise  ratio  (SNR)  are  carried 
out  on  simulated  Raman  spectra  using  the  matched  filter  technique.  The 
matched  filter  is  implemented  using  both  the  DFT  and  the  DHT.  The 
computational  efficiency  and  the  filter  response  in  both  the  methods  are 
compared. 

1.3  Structure   of   the  Thesis 

In  Chapter  2,  we  introduce  the  ET  and  develop  its  properties,  and 
in  Chapter  3  we  define  the  DHT  and  develop  its  properties.  Chapters  4 
and  5  present  a  detailed  development  of  the  radix-2  DIT  and  DIF  algo- 
rithms, respectively.  We  discuss  the  radix-4  DIT  algorithm  in  Chapter 
6,  and  the  split-radix  algorithm  in  Chapter  7.  In  Chapter  8,  we  discuss 
an  application  of  the  DHT  in  the  analysis  of  Raman  spectra,  present  the 
results  of  the  simulation,  and  compare  its  performance  with  the  DFT  in 
that  application.  Finally,  Chapter  9  gives  a  summary  of  the  conclusions 
drawn. 
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2.      THE  HARTLEY  TRANSFORM 

2.1  Introduction 

Since  the  HT  is  so  closely  related  to  the  FT  and  many  of  its 
properties  follow  directly  from  those  of  the  FT,  we  begin  the  chapter  by 
defining  the  FT  of  an  aperiodic  signal.  Then  we  define  the  HT  as 
proposed  by  Bracewell  and  develop  various  properties  of  the  transform. 

2.2  The  Fourier  Transform 

Given  a  signal  x(t),  its  Fourier  transform  X(f)  is  defined  as 


!   (t)e"j2TTft  dt.  (2.1) 


J 


To  recover  x(t)  from  X(f),  we  perform  the  inverse  Fourier  transform 
(IFT)  on  X(f).   The  IFT  is  defined  as 

00 

x(t)  =  |  X(f)ej2"ft  df.  (2.2) 

—00 

Since 

e~j2nft  =  cos(2jift)   -   jsin(2nft), 
we   can  rearrange  Eqn.    (2.1)    as 

00  00 

X(f)   -  J  x(t)cos(2nft)dt  -  j   j  x(t) sin(2nf t)dt .  (2.3) 

— -oo  —00 

As  evident  from  (2.3),  in  general,  the  FT  X(f)  of  a  real  signal  x(t)  is 
complex,  with  the  real  part 
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Re  {X(f)}  =  J  x(t)cos(2nft)dt,  (2.4) 


and  the  imaginary  part 


Im  (X(f)}  =  -   J  x(t)sin(27Tft)dt.  (2.5) 


2.3   The  Hartley  Transform 

The  HT  H(f)  of  a  real  signal  x(t)  is  defined  as 


H(f)  =  J 


x(t)cas(2nft)dt  (2.6) 


where 

cas(2nft)  =  cos(2nft)  +  sin(2nft).  (2.7) 

The  HT  is  symmetric   that  is,  if  we  need  to  recover  x(t)  from  H(f),  we 
perform  the  transformation  in  (2.6)  on  H(f).   Thus  both  the  forward  and 
inverse  transforms  are  identical. 
Using  (2.7)  we  can  rewrite  (2.6)  as 


I 


H(f)  =  I  x(t)cos(2nft)dt  +  I  x(t)sin(2nft)dt.         (2.8) 

_CO  —00 

From    (2.8)    we    observe    that,    unlike    the   FT,    the   HT   of    a    signal    x(t)    is 

real. 

Using  the  trigonometric  identities 

cos(-  x)  =  cos(x) 
and 

sin(-  x)  =  -  sin(x) 
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we  can  compute  II(-f)  as 


H(-f)  =  j  x(t)cos(2nft)dt  -  f  x(t )sin(2nft)dt .         (2.9) 


Now,  by  splitting  H(f)  into  even  and  odd  parts  H  (f)  and  H  (f)  respec- 

c         o 


tively,  we  get 


H  (f)  =  [H(f)  +  H(-f)]  /  2,  (2.10) 

e 


and 


H  (f)  =  [11(f)  -  K(-f)]  /  2.  (2.11) 

o 

Substituting  the  integral  expressions  for  H(f)  and  H(-f)  that  we  ob- 
tained in  (2.8)  and  (2.9)  in  (2.10)  and  (2.11).  we  get 

09 

He(f)  -  J  x(t)cos(2jtft)dt.  (2.12) 


and 


HQ(f)  =  f  x(t)sin(2nft)dt  (2.13) 


2.4      Properties 

The  HT  has  many  interesting  properties  that  can  be  successfully 
employed  in  many  signal  processing  applications.  In  this  section  we 
will   derive    some   of    its    important   properties. 

2.4.1     Relation  Between   the   FT  and   the   HT 

From    (2.4)    and    (2.12),    we   observe    that 
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Re  (X(f)J  =  H  (f). 
e 

and 

from  (2.5)  and  (2.13),  it  is  clear  that 

Im  (X(f)}  =  -  H  (f). 
o 

Therefore,  once  we  have  the  ITT  of  a  signal  x(t),  we  can  compute  its  FT 

as  shown  below. 

X(f)  -  H  (f)  -  j  H  (f).  (2.14) 

e         o 

By  referring  to  (2.4),  (2.5),  (2.12)  and  (2.13)  again,  we  see  that  the 

FT  of  a  signal  z(t)  can  be  computed  from  its  FT  X(f)  as  follows: 

H(f)  =  Re  {X(f)J  -  Im  (X(f)}.  (2.15) 

2.4.2  Linearity 

The  in  is  a  linear  transform. 
Let  us  find  the  transform  of  the  signal  {x  (t)  +  x  (t)}.   Let  the  HT  of 

x  (t)  be  H  (f)  and  that  of  x  (t)  be  H  (f).   Let  the  HT  of 

{xx(t)  +  x2(t))  be  H(f).   Then  n(f)  is  given  by 


H(f)  =  J  [xx(t)  +  x2(t)]cas(2nft)dt. 


The  above  integral  can  be  split  into  two  integrals  as  shown  below 


H(f)  =  J  x1(t)cas(2nft)dt  +  f  x2(t)cas(2nft)dt 


=  H  (£)  +  H2(f). 
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Therefore,  the  HT  of  a  sum  of  signals  equals  the  sum  of  the  HTs  of 
individual  signals. 

2.4.3   Reversal 

If  the  HT  of  x(t)  is  H(f)  then  the  HT  of  x(-t)  is  H(-f ) . 
The  KT  of  x(-t)  is  given  by 


■f-«- 


HT  {x(-t)}  =  J  s(-t)cas(2Trft)dt. 

—00 

Now,    let   -t   =  u.      Then    dt    =   -du.       Making    these    substitutions    in    the 
above    integral,    we   get 


■  J«<« 


HT   U(-t)}    =  J    x(u)cas{2nf(-u)}du 

—00 
00 

=  J    x(u)cas{2nu(-f)}du 

— OD 

=  H(-f). 
Therefore,  if  the  HT  of  x(t)  is  H(f),  the  HT  of  x(-t)  is  H(-f) 

2.4.4  Transforms  of  Even  and  Odd  Functions  ■ 
Let  x(t)  be  an  even  function.   i.e., 
x(t)  =  x(-t). 
The  HT  of  x(t)  is  obtained  as 


H(f)  =  f  x(t)cas(2nft)dt. 
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Since    x(t)     is    an    even    function,     replacing   z(t)    by   x(-t)    in   the    above 
integral,    we    get 


=  1 


H(f)    =  J    x(-t)cas(2nft)dt 

—  €0 

By  the  reversal  theorem,  the  right  hand  side  of  the  above  equation  is 
H(-f).   Therefore, 

H(f)  =  B(-f). 
Therefore,  the  HT  of  an  even  function  is  also  even. 
Now,  let  x(t)  be  an  odd  function    i.e.,  x(t)  ■  -  x(-t). 
If  we  denote  the  HT  of  x(t)  by  H(f), 


-J 


H(f)  =  J  x(t)cas(2nft)dt. 

—  00 

Since    x(t)    =    -    x(-t),    replacing   x(t)    by  -   x(-t)    in   the    above    equation, 
we    get 


•"  I 


H(f)    =  "     J    x(-t)cas(2nft)dt. 

— oo 

By  the  reversal  theorem,  the  right  hand  side  of  the  above  equation  is 
-  H(-f).   Thus,  the  ibove  equation  reduces  to 

H(f)  =  -  B(-f). 
Therefore,  the  HT  of  an  odd  function  is  also  odd. 
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2.4.5  The  Infinite  Integral  Theorem 

The  infinite  integral  of  a  signal  x(t)  can  be  computed  very  easily 
from  its  HT.  This  result  follows  directly  from  the  definition  of  the 
HT.   We  know  that 

00 

H(f)  =  J  x(t)cas(2nft)dt. 

—00 

Substituting  f  =  0  in  the  above  equation,  we  get 

00 

H(0)  =  J  x(t)dt,  (2.16) 

—  00 

since 

cas(O)  =  cos(O)  +  sin(O)  =  1. 

2.4.6  Shift  Theorem 

Now,  let  us  find  the  HT  of  x(t  -  T),  which  is  the  signal  x(t) 
shifted  by  T  units. 
We  know  that  the  HT  of  x(t  -  T)  is  given  by 

00 

HT  {x(t  -  T)}  =  f  x(t  -  T)cas(27Tft)dt.  (2.17) 

—oo 

Let  (t  -  T)  -  h.   Then 

t  =  (T  +  h), 
and 

dt  =  dh. 
Using  the  above  substitutions,  we  rearrange  (2.17)  as 
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-\ 


ITT    U(t   -  T)}    =  j    x(h)cas    {2nf(h  +  T)}dh.  (2.18) 

—  00 

Let   us   expand   cas{2nf(h  +  T)}. 

cas{2irf(h  +  T)}    =  cos    {2nf(h  +  T)}    +   sin{2nf(h  +  T)} 

■  cos(2nfh)cos(2n£T)   -    sin(2nfh)  sin(2itfT) 

+  sin(2nfh)cos(2irfT)    +   cos(2nf  h)  sin(2rtfT)  .       (2.19) 

Collecting    terms    in   the  above   expansion,    we    get 

cas{2nf(h  +  T)}    =  cos(2nfT) {cos(2nfh)    +   sin(2jrfh)} 

+  sin(2nfT) {cos(2nfh)  -  sin(2nfh)}. 
Recognizing  that 

cos(2nfh)    +   sin(2nfh)    =  cas(2nfh)# 
and 

cos(2nfh)    -   sin(2nfh)    =  cos{2?rh(-f ) }    +   sin{2nh(-f)} 

=  cas{2nh(-f)}. 
we   can   rewrite    (2.19)    as 

cas{2nf(h  +  T)}    =  cos(2nfT)cas(2nfh) 

+  sin(2jrfT)cas{2nh(-f)}.  (2.20) 

Substituting  the  result  of  (2.20)  in  (2.18),  we  get 


HT  {x(t  -  T)}  =  f  x(h)[cos(2rtfT)cas(2nfh) 

—00 

+  sin(27TfT)cas{2rrh(-f)}]dh. 
Splitting  the  above  integral  into  two  integrals,  we  get 

09 

HT{x(t  -  T)}    =   cos(2jrfT)    f  x(h)cas(2nfh)dh 
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+  s 


in(2nfT)  f  x(h)cas{2nh(-f )}dh 


=  cos(2jrfT)H(f)  +  sin(27tfT)n(-f ) .  (2.21) 

Note  that  if  x(t)  is  even,  and  thus  H(f)  is  even,  then  the  above  result 
reduces  to 

H(f)[cos(2nfT)  +  sin(2nfT)]  =  cas(2nfT)II(f ) . 
and  if  x(t)  is  odd,  and  thus  H(f)  is  odd,  (2.21)  becomes 
H(f)[cos(2jrfT)  -  sin(27tfT)]  =  cas{2nT(-f )}H(f ) . 

2.4.7  Derivative  Theorem 

From  FT  theory,  we  know  that  if  X(f)  is  the  FT  of  the  signal  x(t), 

then  the  FT  of  its  derivative  is  given  by  j2nf  X(f)   i.e., 

FT  {dx(t)/dt}  =  j2nf  X(f). 

From  (2.14),  we  can  rewrite  j2nf  X(f)  as  follows: 

j2nf  X(f)  =  j2nf  [H  (f)  -  jH  (f)] 

e       o 

=  2nf  [H  (f)  +  jH  (£)]. 
o        e 

where  II  (f)  and  H  (f)  are  the  even  and  odd  parts  of  the  HT  H(f)  of  the 
e         o 

signal  x(t) . 

Since    the   above   expression    is    the    FT    of    the    derivative    of    x(t), 

from   (2.15),    the  HT  of   the   derivative    is   obtained   as 

HT   {dx(t)/dt}    =  Re    [FT    {dx(t)/dt}]   -   Im    [FT    {dx(t)/dt}] 


=  2nfH   (f)   -  2nfH    (f).  (2.22) 

o  e 


Using  the  relations 
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ii  (f)  =  [n(f)  -  ii(-f)]  /  2. 

o 
and 

H    (f)    =    [11(f)    +  H(-f)]    /    2. 

e 

we    rewrite    (2.21)    as 

HT   (dx(t)/dt)    =  -  2nf  H(-f). 
Alternatively,    we    can    use    the    shift    theorem    to    prove    the    derivative 
theorem.      We   know   that    the   derivative   dx(t)/dt    is   defined   as 

dx(t)/dt   =  .    Lim.    n    [x(t   +  h)    -  x(t)]    /    h.  (2.23) 

If  H(f)  is  the  HT  of  x(t),  then  by  the  shift  theorem 

HT  {x(t  +  h)}  =  cos(2nfh)H(f)  -  sin(2nfh)H(-f ) . 

Taking  the  HT  on  both  sides  of  (2.23)  and  applying  the  above  result,  we 

get 

HT   {dx(t)/dt}    -     Liin      _    [cos(2Trfh)H(f)    -   sin(2nf h)H(-f ) 
n >   0 

-   H(f)]    /    h. 

Evaluating  the  above  limit  using  L'Hospital's  rule,  we  get 

HT  {dx(t)/dt}  =  -2!fH(-f). 

2.4.8  HT  of  the  n*   Derivative 

The  FT  of  the  n   derivative  of  a  signal  x(t),  whose  FT  is  X(f),  is 

given  by  (j2nf)nX(f). 
By  (2.14), 

(j2nf)nX(f)  =  j"(2nf)n  [H  (f)  -  jH  (f)], 

e        o 

where  H(f)  is  the  HT  of  x(t). 
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Noting    that 


.n  jnn/2 

J      ■   e 


=   cos(nji/2)    +  jsin(nn/2), 
we    rewrite    the   above   equation  as 


FT    {dnx(t)/dtn}    =    (2nf)n[cos(nn/2)    +  jsin(nn/2)] 

t{H(f)    +  H(-f)}    /   2  -   j{H(f)   -  H(-f)}    /   2], 
(2.15)    tells   us   that    if   we    subtract    the    imaginary   part    of    the    above 

expression    from    its    real    part,    we    get   the   HT  of   the   n        derivative   of 
x(t) .      Therefore, 

HT    {dnx(t)/dtn}    =»    (2nf)n[cos(nn/2)H(f)   -   sin(mr/2)H(-f )  ] . 
However,    when  n   is  odd,    the    first    term   in   the   above   expression  vanishes 
and,    when  n   is   even  the    second   term  vanishes.      Therefore, 

HT   {dnx(t)/dtn}    =  -    (2nf)n   sin(nn/2)H(-f)  (n  odd) 

HT   (dnx(t)/dtn)    =        (2nf)n  cos(nn/2)H(f).  (n  even) 


2.4.9  First-Moment  Theorem 

Given  a  signal  x(t),  its  first  moment  is  given  by 

CO 

F.  M.  =  f  tx(t)dt. 

—OB 

Now,  we  know  that  the  HT  H(f)  of  x(t)  is  given  by 

00 

H(f)  =  f  x(t)cas(2nft)dt. 

—00 
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Taking  the  derivative  of  the  above  equation  with  respect  to  f  yields 

00 

dH(f)/df  =  2n  f  x(t)  t[cos  (2nft)  -  sin  (2nft)]  dt. 

—  00 

Evaluating  the  above  expression  at  f  =  0  gives 


dH 


(0)/df  =  2jt  f  x(t)tdt. 


Therefore,  the  first  moment  can  be  found  as 


*  M#  =  J" 


tx(t)dt  =  [dH(0)/df]  /  2ir. 


2.4.10   Second-Moment  Theorem 

Given  a  signal  x(t),  its  second  moment  is  defined  as 


i.  M.  =  J  t2x(t)dt. 


Differentiating  (2.6)  two  times  with  respect  to  f  yields 


H  (f)  =  4n2  J  t2  x(t)  [-cos(2nft)  -  sin  (2nft)]  dt 


Evaluating  the  above  expression  at  f  =  0  gives 


H  (f)  =  -  4n2  J  t2  x(t)dt 


Therefore,  the  second  moment  is 
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i.  M.  =  J  t2  x(t)dt 


-  H  (0)/4n2. 


2.4.11  Centroid 

Given  a  function  x(t),  its  centroid  is  defined  as 

00  09 

[  J  tx(t)dt]  /  I  x(t)dt  =  H'(0)  /  {2n[H(0)]}. 

—00  —  00 

The  result  follows  directly  from  the  infinite  integral  and  first-moment 
theorems . 

2.4.12  Autocorrelation  Theorem 

From  FT  theory,  we  know  that  the  FT  of  the  autocorrelation  of  x(t) 

is  given  by  |X(f)|   where,  X(f)  is  the  FT  of  x(t). 


However, 


|X(f)|2  =  Re  [X(f)]2  +  Im  [X(f)]2. 


From  (2.14),  it  follows  that 


Since 


and 


|X(f)|2  =  H2(f)  +  H2(f). 
e       o 


H  (f)  =  [H(f)  +  H(-f)]  /  2, 

e 


H  (f)  =  [H(f)  -  H(-f)]  /  2, 
o 
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the  above  equation  can  be  rearranged  as 

|X(f)|2  =  (l/2)[H2(f)  +  H2(-f)]. 
Thus,  the  HT  of  the  autocorrelation  function  is  a  non-negative  and  even 
function.   Observe  that  if  H(f)  is  even,  the  above  expression  becomes 

ET(f),  which  resembles  the  familiar  FT  result. 

2.4.13   Convolution  Theorem 

If  the  FTs  of  signals  x  (t)  and  x  (t)  are  X  (f)  and  X -(f).  respec- 
tively, then  we  know  that  the  FT  of  the  convolution  of  the  two  signals 
is  the  product  of  the  individual  FTs  X.(f)  and  X_(f)   i.e., 

FT  {x^t)  *  x2(t)}  =   X1(f)X2(f) 
By  (2.14),  we  can  rewrite  X  (f)X  (f)  as 

X1(f)X,(f)  =  [H,  (f)  -  jH,  (f)][H_  (f)  -  jEL  (f)l. 

12        le        lo      2e        2o 


where 


and 


H,  (f)  is  the  even  part  of  the  HT  H,  (f)  of  x„(t), 
le  11 

H,  (f)  is  the  odd  part  of  H, (f), 
lo  1 

H.  (f)  is  the  even  part  of  the  HT  H.(f)  of  x.(t), 
Ze  2        2 


H„  (f)  is  the  odd  part  of  H„(f). 
~o  2 


Now,  carrying  out  the  multiplication  in  the  above  equation,  we  get 

X^m.m  =  [H,  (f)H_  (f)  -  H,  (f)H„  (f)] 
12        le    2e       lo    2o 

-  jtH,  (f)H.  (f)  +  H,  (f)H_  (£)]. 
le    2o       lo    2e 


1  n 
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From    (2.15),    the   IIT  of   the   convolution   is   given  by 

HT    {Xl(t)*x,(t)}    =  H,     (f)H.    (f)    -  B,     (f)B.    (f)    +  H      (f)H      (f) 
12  le  2e  lo  2o  le  2o 

+  BL     (f)BL    (f). 
lo  2e 

Collecting  the  terms,  we  get 

HT{x1(t)*x.(t)}    =  H.    (f)    [H.    (f)    +  BL    (f)] 

12  le  Ze  2o 

-  BL    (f)    [H-n(f)    -   B.    (£)]. 

lo  iv  Ze 

=  BL     (f)B0(f)    -  BL     (f)    [B.    (f)    -  B.    (f)]. 
le  2  lo  2o  2e 

=    [{B1(f)    +  B1(-f)}    /   2]B2(f) 

-  [{B1(f)    -   B1(-f)}    /    2]    [{B2(f)    -  B2(-f)} 

-  {B2(f)    +  B2(-f)}]    /    2. 

Simplifying   the   above   expression,    we    get 

BT    {x1(t)*x2(t)}    =    (1/2)    [B1(f)B2(f)    +  B1(-f)B2(f)    +  H    (f)H    (-f) 

-  B1(-f)B2(-f)]. 

Bowever,  if  we  have  any  symmetries  in  the  functions  being  convolved,  the 
above  result  simplifies.   Following  is  a  summary: 

SYMMETRY  BT  OF  TBE  CONVOLUTION 

Bx(f)  is  even  H  <f)H2(f> 

B2(f)  is  even  H  (f)H2(f) 

Both  are  even  B  (f)B  (f) 

B1(f)  is  odd  B1(f)E2(-f) 

B2(f)  is  odd  B1(-f)B2(f) 

Eoth  are  odd  -H  (f)H  (f) 
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2.4.14      Product   Theorem 

If    two    signals    x..(t)    and    x.(t)    have   FTs   X    (f)    and   X -(f).    respec- 
tively,   we   know    that    the   FT   of    the    product    of    the    two    signals     is    the 
convolution  of    the    individual   FTs         i.e., 
FT(x    (t)x2(t)J    =  X1(f)*X2(f). 


By    (2.14), 


X1(f)*X.(f)    =    [H.     (f)    -   jH.     (f)]    •    [H      (f)    -   jH      (f)]. 
12  le  lo  Ze  Zo 


Rearranging    the    right    hand    side    of    the    above    equation    into    real    and   even 

parts,    we    get 

X1(f)*X.(f)    =    [H,    (f)    *  H,    (f)   -  E\    (f)    *   H.    (f)] 
12  le  2e  lo  2o 

-j[H,    (f)    *  H.    (f)    +  H,     (f)    ♦   H,    (f)]. 
le  2o  lo  2e 

(2.15)    tells    us    that    the    HT    of    the   product    is   obtained  by    subtracting 

the    imaginary   part    of    the    above    expression    from    its    real    part. 

Therefore, 

HTU^Ox.U)}    =  H.    (f)    •  H.    (f)    -  H,    (f)    *  H.    (f) 
l  z  le  ze  lo  Zo 

+  h,   (f)  *  n.  (f)  +  n.   (f)  *  n„  (£). 

le  2o  lo  2e 


since 


H.     (f)    =    [H,(f)    +  H\(-f)]    /    2, 
le  1  1 


H.    (f)    =    [H.(f)   -  n\(-f)]    I   2, 
lo  1  1 


H,    (f)    =    [H_(f)    +  H.(-f)]    /   2, 
ze  z  z 


H,    (f)    =    [H.(f)    -   H.(-f)]    /   2, 
Zo  z  z 


the   above   equation  can  be    simplified   as 


HT{x1(t)x2(t)}    =    (l/2)[H1(f)    *   n2(f) 
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+  H^-f)   *  n2(f)   +  Hx(f)   *  H2(-f) 

-  ^(-f)  *  n2(-f)]. 

However,     if   the   functions  being  multiplied  possess    some    symmetries,    the 
above   result    simplifies.      Following    is   a   summary: 

SYMMETRY  HT  OF  THE  PRODUCT 

H^f)    is   even  H   (f)*H    (f) 

H2(f)    is   even  H    (f)*H    (f) 

Both  are   even  H   (f)*H    (f) 

H   (f)    is   odd  H1(f)*H2(-f) 

H2(f)    is  odd  H1(-f)*H2(f) 

Both  are   odd  -H    (f)*E    (f) 


2.4.15      Cross-Correlation  Theorem 

If     two     signals     x     (t)     and     *2(t)     have    FTs    X     (f)     and 

X_(f),    respectively,    the  FT  of   the   cross-correlation  of   the    two    signals 

* 
is    given   by   X    (f)X.(f),    where    the    superscript    '    *    '    indicates   the   com- 
plex  conjugate.      Using   the   relation  between  the  FT  and   the   HT,    we   get 

X1*(f)X.(f)   =    [H.    (f)   -  jH,    (f)]*    [H.    (f)    -   jH.    (f)] 
12  le  lo  2e  2o 

Carrying  out  the  multiplication  and  rearranging  the  above  expression,  we 

get 

X*(f)X.(f)  =  [H,  (f)H.  (f)  +  H,  (f)H.  <f>] 
12        le    2e       lo    2o 
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-   J[Hle(f)H2o(f)    -   Hlo(f)U2e(f)]. 

From    (2.15),    we    know    that    the   IIT  of    the    cross-correlation    is    obtained   by 

subtracting    the     imaginary    part    of    the    above    expression   from    its    real 

part. 

HT{cross-correlation  of   x    (t)    and   x_(t)}    ■ 

H,     (f)H.     (f)    +  H.     (f)H,    (f)    +  H,     (f)H      (f)    -   IT      (f)H„    (f). 
le  2e  lo  2o  le  2o  lo  2e 

Substituting    for    H.     (f),    H,     (f),    H.     (f),    H„     (f)     in    terms    of    IL(f), 

le       lo       2  e       1 o  1 

H.(f),  H  (-f)  and  H_(-f)  in  the  above  expression,  we  get 

HT{cross  correlation  of  x  (t)  and  x  (t)}  = 

(1/2)[H  (£)H  (f)  +  H1(-f)H2(f)  -  H  (f)H2(-f)  +  B^-f )K2(-f )1 . 

However,    if   one   of   the   functions    in  the   cross-correlation  possesses    some 
symmetry,     the    above    result    simplifies.        Following     is     a     summary: 
SYMMETRY  ITT  OF  THE  CROSS   CORRELATION 

H^f)    is   even  H    (f)H    (f) 

H2(f)  is  even  H^-fJH^f) 

Both  are  even  H  (f)B  (f) 

Hx(f)  is  odd  -H  (f)H2(-f) 

E2(f)  is  odd  n1(f)H2(f) 

Both  are  odd  H  (f)H2(f) 


2.4.16  Parseval's  Theorem 

Parseval's  theorem  states  that 
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J  I  x(t)|2dt   =  J  I  X(f)|2df 

—  00  —CO 

where 

X(f)  is  the  FT  of  the  signal  x(t). 

The  above  theorem  holds  true  for  the  HT  also. 

Noting  that 


-I 


x(t)  =   H(f)cas(2nft)df, 


we  can  write 


f  x2(t)dt  =  J"  x(t)  f  H(f)cas(2nft)dfdt. 

■■09  —CO  GO 

Interchanging  the  order  of  integration  on  the  right  hand  side,  we  get 

00  CO  CO 

f  x2(t)dt  =  J  H(f)  f  x(t)cas(2nft)dt  df. 


—  00  —00 


The  second  integral  on  the  right  hand  side  is  the  HT  of  the  signal. 
Therefore, 


f  x2(t)dt  =  j  H2(f)df. 


Wang  [5]  described  a  real  series  representation  for  periodic  signals. 
Since  from  a  signal  processing  point  of  view  we  are  more  interested  in 
data  sequences  than  in  CT  signals,  in  the  next  chapter  the  DHT  is 
defined   and   its   properties   are   derived. 
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3.   THE  DISCRETE  HARTLEY  TRANSFORM 

3.1  Introduction 

We  begin  the  chapter  by  defining  the  DFT  of  a  finite-length 
sequence.   The  DHT  as  defined  by  Bracewell  [2]  is  presented  next. 
Finally,  many  useful  properties  of  the  DIIT  are  derived  and  comparisons 
made  with  those  of  the  DFT. 

3.2  The  Discrete  Fourier  Transform 

The  DFT  {X(k)}  of  a  finite-duration  sequence 

U(n)}  =  {x(0),  x(l),  ...»  x(N  -  1)}  is  defined  as 

N-l 
X(k)  =  (1/N)   ^  x(n)W~nk.  k  =  0.  1,  ...,  N  -  1  (3.1) 

n=0 

where 

W.T  is  called  the  kernel  and  is  given  by  W.,  =  e 
N  N 

Using  the  definition  for  the  kernel  given  above,  we  can  expand  (3.1)  as 

N-l 
X(k)  =  (1/N)  5  x(n)cos(2nnk/N) 
n=0 

N-l 
-  j(l/N)  J  x(n)sin(2;tnk/N),  k  =  0,  1,  ....  N-l.  (3.2) 

n=0 

From  (3.2)  it  is  clear  that  the  DFT  of  a  sequence  is,  in  general, 

complex,  with  the  real  part 
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N-l 
Re  {X(k)}  =  (1/N)  ^  x(n)cos(2jrkn/N).  (3.3) 

n=0 

and  the  imaginary  part 

N-l 
Im  {X(k)}  =  (1/N)  5  x(n)sin(2nkn/N).  (3.4) 

n=0 

To  recover  (x(n)}  from  (X(k)},  we  perform  the  inverse  DFT  (IDFT),  which 

is  defined  as 

N-l 
x(n)  =   ^  X(k)W  nk,  n  =  0,  1,  ...,  N-l  (3.5) 

n=0 

where,  W.t  is  the  kernel  as  defined  earlier. 

N 


3.3   The  Discrete  Hartley  Transform 

The  DHT  H(k)  of  an  N-point  real  sequence  {x(n)}  is  defined  as 

N-l 
H(k)  =  (1/N)  ^  x(n)cas(2nnk/N),  k  =  0.  1,  ....  N-l.  (3.6) 

n=0 

where 

cas(2nnk/N)    =   cos(2nnk/N)    +   sin(2nnk/N). 

Using    the    above    expansion    for    cas (2nnk/N) ,    we    split    (3.6)    into   two 

summations   as   shown  below. 

N-l  N-l 

H(k)    =    (1/N)   ^     x(n)cos(2jrnk/N)    +    (1/N)   ^  x(n)sin(2nnk/N) . 
n=0  n=0 

k  =  0,  1.  ...,  N  -  1.(3.7) 
We  can  compute  H.(-k)  from  the  above  equation  as 
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N-l 
H(-k)  =  (1/N)  ^  x(n)cos{2nn(-k)/N} 
n=0 

N-l 
+  (1/N)  ^  x(n)sin{2nn(-k)/N}.  k  =  0,  1,  ....  N-l.         (3.8) 
n=0 

Using  the  trigonometric  identities 

cos{2nn(-k)}  =  cos(2nnk) 
and 

sin{2nn(-k)}  =  -  sin(2nnk) 

we  can  rewrite  (3.8)  as 

N-l 
H(-k)  =  (1/N)  J  x(n)cos(2nnk/N) 
n=0 

N-l 
-  (1/N)  ^  x(n)sin(2nnk/N),  k  =  0,  1,  ....  N-l.  (3.9) 

n=0 

We  can  split  H(k)  into  even  and  odd  parts,  with  the  even  part  TI  (k) 

B 

being  given  by 

H  (k)  =  [H(k)  +  H(-k)]  /  2,  (3.10) 

and  the  odd  part  H  (k)  obtained  as 
o 

H  (k)  =  [H(k)  -  H(-k)]  /  2.  (3.11) 

o 

Substituting  the  summations  of  (3.7)  and  (3.9)  for  II  ( k )  and  H(-k)  in 

(3.10)  and  (3.11),  we  get 

N-l 
H  (k)  =  (1/N)  ^  x(n)cos(2nnk/N),  k  =  0,  1,  ....  N-l  (3.12) 

n=0 
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N-l 
H  (k)  =  (1/N)  )  x(n)sin(2nnk/N),  k  =  0,  1,  ....  N-l.  (3.13) 

O  L 

n=0 
In  all  the  above  equations,  H(-k)  =  H(N  -  k) .   For  example,  for  a  six- 
teen point  sequence,  i.e.  N  =  16,  H(-l)  is  H(15),  H(-2)  is  H(14)  and  so 
on.   The  notation  H(-k)  instead  of  H(N  -  k)  is  only  an  extension  of  the 
habit  formed  in  dealing  with  the  CT  functions. 

Unlike  the  DFT,  the  DHT  does  not  have  a  separate  inverse  DHT. 
If  we  want  to  recover  the  original  data  sequence  {x(n)},  we  only  need  to 
perform  the  DHT  operation  on  {H(k)J. 

N-l 
x(n)  =  Y   R(t)cas(2nnk/N).  n  =  0,  1,  ....  N-l.  (3.14) 

k=0 

Note    the   absence   of   the    (1/N)    factor   in    the    above    equation   before    the 

summation    symbol,    which  was   present  when  we   defined  H(k).      Its   absence 

does  not   alter   the    fact   that   the  DHT   is    symmetric,    because    by    symmetry 

we    mean    invariance   of   the   kernel    in    (3.6)    and    (3.14).      The   factor    (1/N) 

is   only  a   constant   that   does  not    affect    the   nature   of   the   transform.      In 

fact,    some   authors   omit    (1/N)    in    (3.6)    as  well. 

If  we   compare    (3.7)    and    (3.2),     it    is    evident    that    for    any    real 

sequence    the    DHT    generates    a    real-valued    sequence,    whereas    the  DFT 

produces   a   complex-valued    sequence.       In    fact,     it    is    this    particular 

property  of   the   DHT  that  makes    it   attractive    in  many  applications. 

3.4      Properties 

The  DHT  possesses  many  useful  properties  that  enable  us  to  develop 
fast   algorithms   to   compute    it,    and   also   apply   it    successfully   in  various 
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signal  processing  applications  such  as  digital  filtering,  fast  convolu- 
tion, spectral  analysis,  etc.  In  this  section  many  of  its  useful 
properties  are  developed. 

3.4.1  Relation  Between  the  DFT  and  the  DHT 

The  DFT  and  the  DHT  are  so  closely  related  that  we  can  move  easily 
from  one  transform  to  the  other.  By  looking  at  (3.3)  and  (3.12)  we 
observe  that 

Re{X(k)}  =  n  (k) 

e 

and 

from  (3.4)  and  (3.13)  we  note  that 

Im{X(k)}  =  H  (k). 
o 

Therefore,  if  we  have  the  DHT  {H(k)J  of  an  N-point  real-valued  sequence, 

we  can  find  the  DFT  of  the  sequence  as  shown  below. 

X(k)  =  H  (k)  -  jH  (k).  k  =  0,  1,  ....  N  -  1.         (3.15) 
e        o 

Referring  to  (3.3),  (3.4),  (3.12),  and  (3.13)  we  note  that  if  we 

have  the  DFT  {X(f)}  of  a  real-valued  sequence,  we  can  compute  the  DHT 

(H(k)}  as 

H(k)  =  Re{X(k)}  -  Im{X(k)},  k  =  0,  1.  ....  N  -  1.     (3.16) 

3.4.2  Periodicity  of  the  Kernel 

We  know  that  the  DHT  of  an  N-point  real-valued  sequence  is  given  by 

N-l 
H(k)  =  (1/N)  ^  x(n)cas(2nnk/N),  k  =  0,  1,  ....  N-l. 
n=0 
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In  the   above   equation   the   factor   cas(2nnk/N)    is   called   the   kernel.      Let 
us    focus   on   the   kernel  briefly. 

cas[2nn(k  +  N)/N]    =   cos[2nn(k  +  N)/N]    +   sin[2nn(k  +  N)/N] 
Expanding    the    right    hand    side    using    the    familiar   trigonometric    iden- 
tities,   we    get 
cas[2nn(k  +  N)/N]    =  cos(2jtnk/N)cos(2jinN/N)    -   sin(2nnk/N) sin(2nnN/N) 

+   sin(27tnk/N)cos(2nnN/N)    +   cos(2nnk/N)sin(2jinN/N) . 

-   cos(2jrnk/N)    +   sin(2nnk/N) 

=  cas(2nnk/N) . 

The  above  equation  shows  that  the  kernel  is  periodic  on  N.   In  the  light 

of  this  fact  let  us  compute  the  (k  +  N)th  element  in  the  DHT  of  the  data 

sequence  {x(n) } . 

N-l 
!(k  +  N)  =  (1/N)  J  x(n)cas[2nn(k  +  N)/N],  k  =  0,  1,  ....  N-l. 


HI 

n=0 


Since  the  kernel  is  periodic  on  N, 

N-l 
H(k  +  N)  =  (1/N)  ^  x(n)cas(2jrnk/N),  k  =  0,  1,  ....  N-l. 
n=0 

-  R(k). 

Therefore,  each  of  the  first  N  elements  in  the  DIIT  sequence  repeats 

itself  after  an  interval  of  N  points.   Now,  let  us  compute  the  (n  +  N) 

element  in  the  data  sequence  (x(n)}  from  its  DHT. 

N-l 
x(n  +  N)  =  5  H(k)cas[2rrk(n  +  N)/N],  n  =  0,  1,  ....  N-l. 
k=0 

Since  the  kernel  is  periodic  on  N,  the  above  equation  reduces  to 
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N-l 

x(n  +  N)  =  J  H(k)cas(2nnk/N),  n  =  0,  1,  ....  N-l. 
k=0 

=  x(n). 
As  we  have  seen  above,  the  periodicity  of  the  kernel  imposes  periodicity 
on  the  data  sequence  as  well  as  its  DHT.  If  the  length  of  the  sequence 
is  N,  whenever  we  perforin  the  DHT  on  it,  the  sequence  is  implicitly 
assumed  to  be  periodic  with  period  N.  In  this  respect  the  DHT  is  iden- 
tical to  the  DFT. 

3.4.3   Orthogonality  of  the  Kernel 

The  kernel  of  the  DHT  is  orthogonal  and  real,  which  makes  the  DHT 

an  orthogonal  transform.   Let  us  consider  the  following  summation: 

N-l 

^  cas(2nnk/N)cas(2nmk/N)  (3.17) 

k=0 

Expanding  the  terms  in  the  summation,  we  get 

N-l 

^   [cos(2nnk/N)    +   sin(2nnk/N)l    [cos(2nmk/N)    +   sin(2rrmk/N)  ] 
k=0 

Simplifying   the   above    summation   further,    we    get 

N-l  N-l 

^  cos(27ink/N)cos(2nmk/N)    +     ^   sin(2nnk/N)cos(2mnk/N) 
k=0  k=0 

N-l  N-l 

+     2  cos(2nnk/N)sin(27rmk/N)    +     )   sin(2nnk/N)  sin(2rtmk/N)  . 
k=0  k=0 
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Collecting  the  terms  and  combining  them  using  the  familiar  trigonometric 

identities,  we  get 

N-l  N-l  N-l 

^  cas(2nnk/N)cas(2nmk/N)  =  ^  cos[2n(n  -  m)k/N]  +  ^  sin[2n(n  +  m)k/N]. 
k=0  k=0  k=0 

However,  both  the  summations  on  the  right  hand  side  are  zero. 

Therefore 

N-l 

)  cas(2nmk/N)cas(2nnk/N)  =  0. 
k=0 

Let  m  =  n  in  (3.17).   Then  it  becomes 

N-l  N-l 

^  cas(2Trnk/N)cas(2nnk/N)    =     ^   [cas(2nnk/N)]2 
k=0  k=0 

N-l 


=  ]!  [cos(2nnk/N)  +  sin(2nnk/N) ]' 


k=0 

N-l 
=  ^  cos2(2nnk/N)  +  sin2(2nnk/N) 
k=0 

+  2cos(2nnk/N)sin(2nnk/N) . 

N-l        N-l 
=  J  1  +2   ^  cos(2nnk/N)sin(2nnk/N) 
k=0        k=0 

=  N  +0  (Since  sine  and  cosine  or  orthogonal 

to  each  other  over  one  full  cycle) 

=  N. 

Therefore,  we  have  proved  that 
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N-l 

)  cas(2nmk/N)cas(27tnk/N)    =0,     if   m    is   not    equal    to   n   and 
k=0 

=  N,     if   m  =   n. 

Hence,  the  kernel  is  orthogonal  and  the  DHT  is  an  orthogonal  transform. 

The  DFT  also  is  an  orthogonal  transform  but  the  kernel  is  complex  in  its 

case . 

3.4.4   Linearity 

The  DHT  is  a  linear  transform  like  the  DFT.   Let  us  consider  two 
N-point  sequences  {x  (n)}  and  {x_(n)}.   Then  the  DHT  of  the  summation  is 

given  by 

N-l 
H(k)  -  (1/N)   ^   [x1(n)  +  x2(n)]cas(2nnk/N),  k  =  0.  1.  ....  N-l. 
n=0 

However,  we  can  split  the  above  summation  into  two  summations  as  shown 

N-l  N-l 

H(k)  =  (1/N)   ^  x1(n)cas(2jrnk/N)  +  (1/N)   ^  x2  (n)cas(2nnk/N) , 
n=0  n=0 

k  =  0.  1,  ....  N  -  1. 

=  Hx(k)  +  H2(k) 

Where,  H  (k)  and  H  (k)  are  the  DHTs  of  the  sequences  {x  (n)}  and 

{x  (n)}  respectively.   Therefore,  the  DHT  of  a  sum  of  sequences  is  the 

the  sum  of  the  DHTs  of  the  sequences  taken  independently.   If  lengths  of 
the  sequences  are  not  identical,  we  make  them  so  by  adding  zeros  to  the 
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data  sequence  having  less  number  of  data  points  before  applying  this 
theorem. 

3.4.5  Reversal 

Let  us  take  an  N-point  data  sequence  {x(n)J  and  arrange  all  the 
elements  in  the  reverse  order  to  get  the  sequence  {x(N  -  n)}.  The  DHT 
of  {x(N  -  )}  is  simply  a  reversed  version  of  the  DHT  {E(k)}  of  {x(n)}, 
except  for  the  first  element.  That  is,  if  {x(0),  x(l),  ....  x(7)}  has 
the  DIIT  sequence  {H(0),  H(l),  ....  H(7)},  then  the  DHT  of  the  sequence 
U(0),  x(7),  x(6),  ....  x(l)}  is  {H(0),  H<7).  H(6),  ...,  H(l)}. 

3.4.6  Transforms  of  Even  and  Odd  Sequences 

We  know  that  the  DHT  of  a  sequence  {x(n)}  is  given  by 

N-l 
H(k)  =  (1/N)   5  x(n)cas(2nnk/N),  k  =  0,  1,  ....  N-l. 
n=0 

For    a    data    sequence    that    is    even,    we   can   replace   x(n)   by  x(-n)    in  the 

above   equation. 

N-l 
H(k)    =    (1/N)      ^     x(-n)cas(2nnk/N),    k  =   0,    1,    ...,    N-l. 
n=0 

By  the   reversal    theorem,    the   right   hand    side    of    the    above    equation    is 

H(-k). 

Therefore 

H(k)    =  H(-k). 

The  above  result  shows  that  the  DHT  of  an  even  real-valued  sequence 

is  also  even. 


The  above  result  can  also  be  proved  via  the  DFT.   From  (3.15)  we 

know  that  X(k)  is  related  to  H(k)  by 

X(k)  =  H  (k)  -  jH  (k). 
e        o 

Since  the  DFT  of  a  real  and  even-sequence  is  real  and,  as  such,  does  not 

have  any  imaginary  part,  from  the  above  equation  we  conclude  that  the 

DIIT  of  the  data  sequence  does  not  have  any  odd  part. 

Now,  let  us  consider  a  data  sequence  that  is  odd.   i.e., 

x(n)  =  -  x(-n) . 

The  DHT  of  {x(n)}  is  given  by, 

N-l 
H(k)  =  (1/N)   ^  x(n)cas(2nnk/N),  k  =  0,  1,  ...,  N-l. 
n=0 

Since  the  data  sequence  is  odd,  we  replace  x(n)  by  -x(-n)  in  the  above 

equation. 

N-l 
H(k)  =  (1/N)   ^  -x(-n)cas(2jrnk/N). 
n=0 

By  the  reversal  theorem,  the  right  hand  side  is  -H(-k) .   Hence 

H(k)  =  -  H(-k). 

We    can    reach    the    same    conclusion   from  our   knowledge    that    the   DFT   of   a 

real    and  odd-sequence    is    imaginary.      Applying   that    fact    to  Eqn.     (3.15), 

we    can  conclude    that    the  DHT  of   an  odd   and   real    sequence    is   also   odd. 

3.4.7      Shift  Theorem 

Let    us    consider    a   data    sequence    (x(n)}    whose   DHT    is    (H(k)}.      Now, 
let   us    circularly    shift    the    data    sequence    by   m    samples    to    the    left. 
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When  we  do  that,  the  samples  that  are  pushed  beyond  the  left  edge  of  the 
sequence  will  reappear  at  its  right  end,  thus  always  maintaining  the 
length  of  the  data  sequence  as  N. 
For  example,  consider  an  8-point  data  sequence  {x(n)}  as  shown  below. 

(x(0),  x(l),  x(2),  x(3),  x(4),  x(5).  x(6),  x(7)}. 
If  we  circularly  shift  {x(n)}  by  two  samples  to  the  left,  we  get  the 
following  sequence. 

(x(2),  x(3),  x(4),  x(5),  x(6),  x(7)»  x(0),  x(l)}. 

Observe  that  the  length  of  the  data  sequence  is  same,  that  is  eight,  in 

both  the  cases. 

Now,  let  us  find  the  DHT  of  such  a  shifted  sequence. 

N-l 
H[x(n  +  m)]  =  (1/N)   ^  x(n  +  m)cas (2nnk/N) ,  k  =  0,  1,  ...,  N  -  1.(3.18) 

n=0 

where  m  is  the  amount  of  shift. 

Let   n  +  m  =  t.   Then  n  =  t  -  m. 

Also,    when  n  =  0,    t  =  m  and  when  n=N-l,    t  =  m  +  N  -   1. 

Making    the   above    substitutions    for    n    and    the    limits    of    summation    in 

(3.18),    we    get 

N-l+m 
H[x(n  +  m)]   =    (1/N)        ^       x(t)cas[2n(t  -  m)k/N],    k  =  0,    1,    ....    N-l. 

t=m 

Expanding    cas[2n(t   -  m)k/N]    and    splitting   the   above    summation    into   two, 
we   get 

N-l+m 


[x(n  +  m)]    =    (1/N)        ^       x(t)cos{2n(t  -  m)k/N} 


H 

t  =  B! 
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+    (1/N)        ^     x(t)sin(2ji(t   -  m)k/N}»    k  =   0.    1 N  -   1. 

t=m 

Using    the    familiar    trigonometric    identities   we    can    rewrite    the    above 

summations   as    follows: 

N-l+m 
H[x(n   +  m)]    =    (1/N)        )     x(t) [cos(2ntk/N)cos (2nmk/N) 

t=m 

+  sin(2ntk/N)sin(2nmk/N)] 

N-l+m 
+  (1/N)    ^  x(t)[sin(2ntk/N)cos(2nmk/N) 
t=m 

-  cos(2ntk/N)sin(2nmk/N)]. 

Collecting  the  terms,  we  get 

N-m+1 
H[x(n  +  m)]  =  (l/N)cos(2mnk/N)    5  x(t) [co$(2ntk/N)  +  sin(2ntk/N) ] 

t=m 

N-m+1 
-  (1/N)sin(27rmk/N)    ^  x(t )  [cos(2ntk/N)  -  sin(2ntk/N)  ]  . 

t=m 

Recognizing  that 

cos(2ntk/N)  +  sin(2ntk/N)  =  cas(2ntk/N) 

cos(2ntk/N)  -  sin(2ntk/N)  =  cas{2nr(-k)/N) , 

we  can  rewrite  the  above  equation  as 

N-l+m 
H[x(n  +  m)]  =  cos(2jtmk/N)(l/N)    ^  x(t )cas(2ntk/N) 

t=m 
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N-l+m 
-  sin(2jrmk/N)(l/N)    J  x(t)cas{2nt (-k)/NJ . 

t=m 

H[x(n  +  m)]  =  cos(2nmk/N)H(k)  -  sin(2jrmk/N)H(-k) 

=  cos(2nmk/N)H(k)  -  sin(2jrmk/N)II(N-k) . 

If  H(k)  is  even,  then  the  above  result  becomes 

H(k)[cas{2nm(-k)}], 

and  if  H(k)  is  odd,  it  will  be  H(k)cas(2mnk/N) . 

3.4.8   Convolution 

Circular  convolution  of  two  data  sequences  can  be  implemented 
easily  using  the  DHT.  To  perform  circular  convolution  of  two  sequences, 
both  the  sequences  must  be  of  identical  length,  unlike  in  linear  con- 
volution where  differing  lengths  are  allowed.  Circular  convolution 
involves  time  reversing  one  sequence,  overlaying  it  on  the  other  and 
then  carrying  out  the  convolution  in  the  usual  manner,  remembering  that 
whenever  we  shift  an  element,  the  shift  is  circular. 

If  we  have  two  data  sequences  {x  (n)}  and  {x-(n)},  each  of  length 

N,  then  their  circular  convolution  is 

N-l 
y(n)  =  (1/N)    J  x1(h)x2(n  -  h) ,  n  =  0.  1,  ...,  N-l. 
h=0 

If  the  DFTs  of  {x  (n)}  and  {x  (n)}  are  X  (k)  and  X  (k)  respec- 
tively, then  the  DFT  of  the  circular  convolution  of  the  two  sequences  is 
F[x1(n)  *  *2(n)]  =  X^kJX^k). 

By  Eqn.  (3.15),  we  rewrite  the  right  hand  side  as 
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X,  (k)X.(k)  =  [n,   (k)  -  jil   (k)][n_   (k)  -  jn_   (k)] 

l  z  le  lo  it  zo 

Where    E,     (k),    B„     (k),    fl„     (k),    B„     (k)    are    the    even-    and    odd-    parts 
le  lo  2e  2o 

respectively,    of   the  DHTs   B    (k)    and  H    (k)    of   the   data    sequences. 

Carrying   out    the   multiplication    in   the    above    equation,    we    get 

X,(k)X„(k)    =    [B,     (k)H0    (k)    -   B,     (k)E,    (k)] 
12  le  2e  lo  2o 

-j[E,    (k)B.    (k)    +  H,    (k)E.    (k)]. 
le  2o  lo  2e 

Using   the   result   of    (3.16)   we   obtain   the   DUT  of   the    convolution  as 

E[x,(n)    •   x.(n)]    =  B.    (k)[H.    (k)    +  B„    (k)]    -  B,    (k)[B0    (k)    -  B,    (k)]. 
12  le  2e  2o  lo  2o  2e 

Writing  B,  (k),  B.  (k),  B.  (k),  B.  (k)  in  terms  of  Hfk),  Bf-k),  fl„(k), 
le      2e      lo      2o  112 

and  B  (-k) ,  we  get 

B[x1(n)*x2(n)]  =  (1/2)  [B1(k)B2  (k)  +  H  <-k)H  <k)  +  B1(k)B2(-k) 

-  B1(-k)B2(-k)]. 

In  some  special  cases  we  get  simpler  results.   Following  is  a  summary: 
SYMMETRY  PET  OF  TEE  CIRCULAR  CONVOLUTION 

(x^n)}  is  even  H  (k)H  (k) 

(x2(n)}  is  even  H  (k)H  (k) 

Both  sequences  are  even  B  (k)fl.(k) 

(x1(n))  is  odd  B  (k)B2(-k) 

U2(n)}  is  odd  B1(-k)E2(k) 

Both  sequences  are  odd  -  B  (k)E  (k) 


V/henever  we  need  to  perform  linear  convolution  of  two  sequences,  {x  (n)} 

of  length  N  and  {x  (n)}  of  length  N  ,  N  >  N  .  we  add  a  string  of  zeros 

to  both  the  sequences  until  the  length  of  each  becomes  at  least  equal  to 
(N  +  N.  -  1),  and  then  perform  circular  convolution  on  the  resulting 

sequences. 


3.4.9      Product  Theorem 

Let    us    consider   two   data    sequences    {x    (n)}     and    (x.(a)},     each    of 

length   N.       If   their  DFT   sequences    are    {X    (k)}    and    {X    (k)}    respectively, 

then   the   DFT  of   their  product    is    given  by    the    circular    convolution    of 
X    (k)    and  X    (k) .   That    is, 

F[x1(n)x2(n)]    =    (1/N)    [X^k)    *  X2(k)] 

Using    the    result    of    (3.15)    we    can    rewrite    the   r *. _,   t   hand   side    of   the 

above    equation   as 

(l/NmCk)    *   X.(k)]    =    (1/N)[H.     (k)    -   jH,     (k)]    *    [H,    (k)    -   jH_    (k)] 
12  le  lo  2e  lo 

Carrying  out  the  multiplication  and  rearranging  the  above  expression,  we 

get 

(l/NHX^k)  *  X.(k)]  =  (1/N)[{H.  (k)  *  H.  (k)  -  H,  (k)  *  H.  (k)} 
12  le       Ze        lo       lo 

-j{H,  (k)  *  H„  (k)  +  H,  (k)  *  H.  (k)}]. 
le       2o       lo       2e 

Using  the  result  of  (3.16)  we  obtain  the  DHT  of  the  product 

{x  (n)x.(n)}  as  shown  below. 

H[x, (n)x„(n)]  =  (1/N) [H,  (k)*{H.  (k)  +  H.  (k)} 
12  le      Ze       zo 
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h,  (k)*{n.  (k)  -  n  (k)}]. 

lo       2o        2e 


Observing  that 


H,  (k) 
le 


[B ,  (k)  +  IL(-k)]/2   n,  (k) 
11  lo 


[H^k)  -  K  (-k)]/2. 


H2e(k) 


[H,(k)  +  n,(-k)]/2   n.  (k) 


[n2(k)  -  H2(-k)]/2 


we  get 

H[Xl(n)x2(n)]  =  (l/2N)[H1(k)  *  H2(k)  +  Il^-k)  ♦  ^(k) 

+  H  (k)  *  H  <-k)  -  D  (-k)  *  D  (-k)] 

If  there  is  any  symmetry  in  one  or  both  of  the  multiplying  sequences, 
the  above  result  simplifies.   Following  is  a  summary. 


SYMMETRY 


DHT  OF  TITE  PRODUCT 


{x  (n)  is  even} 


{x_ (n)  is  even} 


Both  sequences  are  even 


{x  (n)  is  odd} 


(x.(n)  is  odd} 


Both  sequences  are  odd 


1/N)[H    (k)H2(k)] 


1/N)[H    (k)H    (k)] 


l/N)[H1(k)H2(k)] 


l/N)[H1(k)H2(-k)] 


l/N)[ni(-k)H2(k)} 


l/N)[-H1(k)H2(k)] 


3.4.10     Cross-Correlation  Theorem 

Circular    cross-correlation    of    two    data    sequences     {x    (n)}     and 

{x    (n)}    is   given  by 

N-l 
y(n)    =    (1/N)        ^     x1(h)x2(n  +   h) ,    n  =   0,    1,    ....    N  -   1.  (3.19) 

h=0 
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Collecting  the  terms  and  combining  them  using  the  familiar  trigonometric 

identities,  we  get 

N-l  N-l  N-l 

2  cas(2nnk/N)cas(2nmk/N)  =  ^  cos[2n(n  -  m)k/N]  +  1  sin[2n(n  +  m)k/N]. 
k=0  k=0  k=0 

However,  both  the  summations  on  the  right  hand  side  are  zero. 

Therefore 

N-l 

)  cas(2nmk/N)cas(2nnk/N)  =  0. 
k=0 

Let  m  =  n  in  (3.17).   Then  it  becomes 

N-l  N-l 

2  cas(2nnk/N)cas(2nnk/N)    =     ^    [cas(2nnk/N) ]2 
k=0  k=0 

N-l 
=     ^   [cos(2nnk/N)    +   sin(2jmk/N)  ]2 
k=0 

N-l 

r  2  2 

=      )  cos    (2nnk/N)    +   sin    (2nnk/N) 

k=0 

+  2cos(2nnk/N)sin(2nnk/N). 

N-l        N-l 
=  ]  1  +2   ^  cos(2nnk/N)sin(2nnk/N) 
k=0        k=0 

=  N  +0  (Since  sine  and  cosine  or  orthogonal 

to  each  other  over  one  full  cycle) 

=  N. 

Therefore,  we  have  proved  that 
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M-l 

)  cas(2nmk/N)cas (2nnk/N)  =0,  if  m  is  not  equal  to  n  and 
k=0 

=  N,  if  m  =  n. 

Hence,  the  kernel  is  orthogonal  and  the  DHT  is  an  orthogonal  transform. 

The  DFT  also  is  an  orthogonal  transform  but  the  kernel  is  complex  in  its 

case . 

3.4.4   Linearity 

The  DHT  is  a  linear  transform  like  the  DFT.   Let  us  consider  two 
N-point  sequences  {x1(n)}  and  (x_(n)}.   Then  the  DHT  of  the  summation  is 

given  by 

N-l 
H(k)  =  (1/N)   ])   [x1(n)  +  x2(n)]cas(2nnk/N),  k  =  0,  1.  ...»  N  -  1. 
n=0 

However,  we  can  split  the  above  summation  into  two  summations  as  shown 

N-l  N-l 

H(k)  =  (1/N)   ^  x1(n)cas(2nnk/N)  +  (1/N)   ]>  x2(n)cas(2nnk/N) . 
n=0  n=0 

k  =  0,  1 N-l. 

=  Hx(k)  +  H2(k) 

Where,  H  (k)  and  H  (k)  are  the  DHTs  of  the  sequences  {x  (n)}  and 

(x_(n)}  respectively.   Therefore,  the  DHT  of  a  sum  of  sequences  is  the 

the  sum  of  the  DHTs  of  the  sequences  taken  independently.   If  lengths  of 
the  sequences  are  not  identical,  we  make  them  so  by  adding  zeros  to  the 
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data  sequence  having  less  number  of  data  points  before  applying  this 
theorem. 

3.4.5  Reversal 

Let  us  take  an  N-point  data  sequence  (x(n)}  and  arrange  all  the 
elements  in  the  reverse  order  to  get  the  sequence  {x(N  -  n)}.  The  DHT 
of  {x(N  -  )}  is  simply  a  reversed  version  of  the  DHT  {E(k)}  of  {x(n)}, 
except  for  the  first  element.  That  is,  if  {x(0),  x(l),  ....  x(7)}  has 
the  DHT  sequence  {H(0),  H(l),  ....  H(7)},  then  the  DHT  of  the  sequence 
{x(0),  x(7),  x(6).  ....  x(l)}  is  {H(0),  H(7),  H(6).  ....  H(l)}. 

3.4.6  Transforms  of  Even  and  Odd  Sequences 

We  know  that  the  DHT  of  a  sequence  {x(n)}  is  given  by 

N-l 
H(k)  =  (1/N)   J  x(n)cas(2nnk/N),  k  =  0,  1,  ...,  N-l. 
n=0 

For    a    data    sequence    that    is    even,    we   can  replace   x(n)   by  x(-n)    in  the 

above   equation. 

N-l 
H(k)    =    (1/N)      ^     x(-n)cas(2nnk/N),    k  =   0,    1,    ....    N-l. 
n=0 

By  the   reversal    theorem,    the   right   hand    side    of    the    above    equation    is 

H(-k). 

Therefore 

H(k)    =  H(-k). 

The  above  result  shows  that  the  DHT  of  an  even  real-valued  sequence 

is  also  even. 
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The  above  result  can  also  be  proved  via  the  DFT.   From  (3.15)  we 

know  that  X(k)  is  related  to  H(k)  by 

X(k)  =  H  (k)  -  jH  (k). 
e        o 

Since  the  DFT  of  a  real  and  even-sequence  is  real  and,  as  such,  does  not 

have  any  imaginary  part,  from  the  above  equation  we  conclude  that  the 

DHT  of  the  data  sequence  does  not  have  any  odd  part. 

Now,  let  us  consider  a  data  sequence  that  is  odd.   i.e., 

x(n)  =  -  x(-n) . 

The  DHT  of  (x(n)}  is  given  by, 

N-l 
H(k)  =  (1/N)   ^  x(n)cas(2nnk/N),  k  =  0,  1,  ....  N-l. 
n=0 

Since  the  data  sequence  is  odd,  we  replace  x(n)  by  -x(-n)  in  the  above 

equation. 

N-l 
H(k)  =  (1/N)   ^  ~x(-n)cas(2nnk/N). 
n=0 

By  the  reversal  theorem,  the  right  hand  side  is  -n(-k).   Hence 

H(k)  =  -  H(-k). 

We  can  reach  the  same  conclusion  from  our  knowledge  that  the  DFT  of  a 

real  and  odd-sequence  is  imaginary.   Applying  that  fact  to  Eqn.  (3.15), 

we  can  conclude  that  the  DHT  of  an  odd  and  real  sequence  is  also  odd. 

3.4.7   Shift  Theorem 

Let  us  consider  a  data  sequence  {x(n)}  whose  DHT  is  {H(k)}.   Now, 
let  us  circularly  shift  the  data  sequence  by  m  samples  to  the  left. 


-34- 


When  we  do  that,  the  samples  that  are  pushed  beyond  the  left  edge  of  the 
sequence  will  reappear  at  its  right  end,  thus  always  maintaining  the 
length  of  the  data  sequence  as  N. 
For  example,  consider  an  8-point  data  sequence  {x(n)J  as  shown  below. 

(x(0),  x(l),  x(2),  x(3),  x(4),  x(5),  x(6),  x(7)}. 
If  we  circularly  shift  {x(n)}  by  two  samples  to  the  left,  we  get  the 
following  sequence. 

(x(2),  x(3),  x(4),  x(5),  x(6),  x(7),  x(0),  x(l)}. 

Observe  that  the  length  of  the  data  sequence  is  same,  that  is  eight,  in 

both  the  cases. 

Now,  let  us  find  the  DHT  of  such  a  shifted  sequence. 

N-l 
H[x(n  +  m)]  =  (1/N)   J  x(n  +  «a)cas (2nnk/N) .  k  =  0,  1.  ...,  N  -  1.(3.18) 

n=0 

where  m  is  the  amount  of  shift. 

Let   n  +  m  =  t.   Then  n  =  t  -  m. 

Also,   when  n  =  0,    t  ■  m  and  when  n  =  N  -  1,    t=m+N-l. 

Making   the  above    substitutions    for    n    and    the    limits    of    summation    in 

(3.18),   we  get 

N-l+m 
H[x(n  +  m)]   =    (1/N)        ^       x(t)cas[2n(t  -  m)k/Nj,    k  =  0,    1,    ....    N-l. 

t=m 

Expanding    cas[2n(t  -  m)k/N]    and   splitting   the   above    summation   into  two, 

we   get 

N-l+m 


[x(n  +  m)]   =    (1/N)        ^       x(t)cos{2n(t  -  m)k/N} 


H. 
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+    (1/N)        5      x(t)sin{2n(t   -  m)k/N},    k  =   0,    1,     ....    N  -    1. 
t=m 

Using    the    familiar    trigonometric    identities   we    can    rewrite    the    above 

summations   as   follows: 

N-l+m 
H[x(n  +  m)]    =    (1/N)        )     x(t ) [cos(2ntk/N)cos(2mnk/N) 

t=m 

+  sin(2ntk/N)sin(2nmk/N)] 

N-l+m 
+  (1/N)    2  x(t)[sin(2ntk/N)cos(2nmk/N) 
t=m 

-   cos(2jTtk/N)sin(2nmk/N)]. 

Collecting   the   terms,    we    get 

N-m+1 
H[x(n  +  m)]    =    (l/N)cos(2nmk/N)         >     x(t>  [cos(2irtk/N)    +   sin(2ntk/N)] 

t=m 

N-m+1 
-  (1/N)sin(27tmk/N)    5  x(t ) [cos(2ntk/N)  -  sin(2ntk/N)  ]  . 

t=m 

Recognizing  that 

cos(2ntk/N)  +  sin(2ntk/N)  =  cas(2ntk/N) 

cos(2ntk/N)  -  sin(2ntk/N)  =  cas{2nT(-k)/N} , 

we  can  rewrite  the  above  equation  as 

N-l+m 
H[x(n  +  m)]  =  cos(2nmk/N)(l/N)    J  x(t)cas(2*tk/N) 

t=m 
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N-l+m 
-  sin(27tmk/N)(l/N)    5  x(t )cas{2nt (-k) /N} . 

t=m 

H[x(n  +  m)]  =  cos(2nmk/N)H(k)  -  sin(2nmk/N)H(-k) 

=  cos(2nmk/N)H(k)  -  sin(2nmk/N)n(N-k) . 

If  H(k)  is  even,  then  the  above  result  becomes 

H(k)[cas{2nm(-k)}], 

and  if  H(k)  is  odd,  it  will  be  H(k)cas(2mnk/N) . 

3.4.8   Convolution 

Circular  convolution  of  two  data  sequences  can  be  implemented 
easily  using  the  DHT.  To  perform  circular  convolution  of  two  sequences, 
both  the  sequences  must  be  of  identical  length,  unlike  in  linear  con- 
volution where  differing  lengths  are  allowed.  Circular  convolution 
involves  time  reversing  one  sequence,  overlaying  it  on  the  other  and 
then  carrying  out  the  convolution  in  the  usual  manner,  remembering  that 
whenever  we  shift  an  element,  the  shift  is  circular. 

If  we  have  two  data  sequences  {x  (n)}  and  {z. (n)},  each  of  length 

N,  then  their  circular  convolution  is 

N-l 
y(n)  =  (1/N)    ^  x1(h)x2(n  -  h) ,  n  =  0,  1,  ....  N-l. 
h=0 

If  the  DFTs  of  {x  (n)}  and  {x  (n)}  are  X  (k)  and  X  (k)  respec- 
tively, then  the  DFT  of  the  circular  convolution  of  the  two  sequences  is 
F[x1(n)  *  *2(n)]  =  X^kJX^k). 

By  Eqn.  (3.15),  we  rewrite  the  right  hand  side  as 
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x.ujx.U)  =  [ir   (k)  -  jir   (k)][n,   (k)  -  jn.   (k)] 

12  le  lo  ze  zo 

Where    E,     (k),    H„     (k),    H„     (k).    H„     (k)    are    the    even-    and    odd-    parts 
le  lo  2e  2o 

respectively,    of   the  DHTs  H    (k)    and  H    (k)    of    the   data    sequences. 

Carrying   out    the   multiplication    in   the    above    equation,    we    get 

Xt(k)X,(k)    =    [H.    (k)B.    (k)    -  H.    (k)BL    (k)] 

12  le  Ze  lo  2o 

-jfH,     (k)H,,    (k)    +   H,    (k)H„    (k)]. 
le  2o  lo  2e 

Using    the    result    of    (3.16)    we    obtain    the    DliT   of    the    convolution   as 

H[x,(n)    *   x.(n)]    =  H,    (k) [H„    (k)    +  H.    (k)]    -  E\    (k) [H_    (k)    -  H„    (k)]. 
1  2  le  Ze  2o  lo  2o  2e 

Writing  H.  (k),  H.  (k),  H.  (k),  H,  (k)  in  terms  of  H,(k),  B(-k).  H„ (k), 
le      2e      lo      2o  11       2 

and  H  (-k) ,  we  get 

H[x1(n)*x2(n)]  =  (1/2)  [H  (k)H  (k)  +  H1(-k)H2(k)  +  H  <k)H  <-k) 

-  H  (-k)H  (-k)]. 

In  some  special  cases  we  get  simpler  results.   Following  is  a  summary: 
SYMMETRY  PITT  OF  THE  CIRCULAR  CONVOLUTION 

(x1(n)}  is  even  H  (k)H  <k) 

(x2(n)}  is  even  H  (k)H2(k) 

Both  sequences  are  even  H  (k)H  (k) 

(x1(n))  is  odd  H1(k)H2(-k) 

(x2(n)}  is  odd  H1(-k)H2(k) 

Both  sequences  are  odd  -  H  (k)E  (k) 


Whenever  we  need  to  perform  linear  convolution  of  two  sequences,  {x  (n)} 

of  length  N   and  {x_(n)}  of  length  N  ,  N   >  N  ,  we  add  a  string  of  zeros 

to  both  the  sequences  until  the  length  of  each  becomes  at  least  equal  to 
(N  +  N   -  1),  and  then  perform  circular  convolution  on  the  resulting 

sequences. 


3.4.9      Product  Theorem 

Let   us   consider   two   data    sequences    {x    (n)}    and    {x_(n)},    each   of 

length   N.       If   their  DFT   sequences   are    {X    (k)}    and    {X.(k)J    respectively, 

then  the   DFT  of   their  product    is    given  by    the    circular    convolution    of 
X   (k)    and  X    (k) .   That    is, 

F[x1(n)x2(n)]    =    (1/N)    [X^k)    *  X2(k)] 

Using    the    result    of    (3.15)    we    can    rewrite    the   r^^   t   hand   side   of   the 

above   equation   as 

<1/N)[X, (k)    *  X„(k)]    =    (1/N)[H,     (k)    -   jH,     (k)]    •    [H.     (k)    -   jH.    (k)] 
12  le  lo  Ze  zo 

Carrying  out  the  multiplication  and  rearranging  the  above  expression,  we 

get 

(l/NHX^k)  *  X.(k)]  =  (1/N)[{H.  (k)  *  H_  (k)  -  H,  (k)  *  II.  (k)} 
12  le       Ze       lo       lo 

-j{R\,  (k)  *  H„  (k)  +  H,  (k)  *  H.  (k)}]. 
le       2o       lo       2e 

Using  the  result  of  (3.16)  we  obtain  the  DHT  of  the  product 

{x  (n)x.(n)}  as  shown  below. 

Hfx^nJx.U)]  =  (1/N>[H,  (k)*{H.  (k)  +H   (k)} 
12  le      2e       20 
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-  h,  (k)*{H.  (k)  -  n_  (k)}]. 

lo      2o       Ze 


Observing  that 


H,  (k)  =  [^(k)  +  H,  (-k)]/2  H,  (k)  =  [H(k)  -  ir(-k)]/2, 

le        1       1  lo        1       1 

h„  (k)  =  [H„(k)  +  iL(-k)]/2  n.  (k)  =  [n„(k)  -  H„(-k)]/2 

Ze        Z       Z  Zo        Z       z 


we  get 


H[x  (n)x2(n)]  =  <1/2N)[H  <k>  *  H2(k)  +  ^(-k)  *  O^k) 

+  e  (k)  *  h  <-k)  -  n  (-k)  *  n2(-k)] 

If  there  is  any  symmetry  in  one  or  both  of  the  multiplying  sequences, 

the  above  result  simplifies.   Following  is  a  summary. 

SYMMETRY  DHT  OP  THE  PRODUCT 

{x1(n)  is  even}  (1/N)  [H^UH^k)  ] 

{x2(n)  is  even}  (1/WIH  (k)^(k)] 

Both  sequences  are  even  (1/N) [H  (k)H  (k) ] 

{x1(n)  is  odd}  (l/N)[H1(k)H2(-k)] 

{x2(n)  is  odd}  (1/N)[H  (-k)H2(k)] 

Both   sequences   are   odd  (1/N)[-H    (k)H    (k)] 


3.4.10      Cross-Correlation  Theorem 

Circular    cross-correlation    of    two    data    sequences    {x    (n)}     and 

{x2(n)}    is   given  by 

N-l 

y(n)    =    (1/N)        ^     x^M^*  +   h)'    n  =   °'    lf    •*•'    N  "   l'  (3.19) 

h=0 
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Note  the  similarity  between  the  circular  cross-correlation  and  the 
circular  convolution  of  two  data  sequences.   The  difference  is  that  in 
the  convolution  operation,  we  flip  one  of  the  data  sequences  before 
carrying  out  the  operation,  whereas  in  cross-correlation  we  do  not. 
If  the  DFT  sequences  of  {x  (n)}  and  {x  (n)}  are  {X  (k)}  and  {X  (k)} 

respectively,  then  the  DFT  of  their  cross-correlation  is 

F[Cross  correlation  of  x  (n)  and  x.(n)]  =  X  (k)X  (k) 

* 
Where  X1 (k)  is  the  complex  conjugate  of  X1 (k) . 

* 
Using  the  result  of  (3.15),  we  rewrite  X  (k)X.(k)  in  terms  of  the  DUT  of 

the  cross  correlation  H(k) . 

X*(k)X.(k)  =  [H,  (k)  +  jH.  (k)][H.  (k)  -  jE„  (k)] 
12        le        lo      2e        2o 

Where  H,  (k),  H,  (k),  H„  (k),  H„  (k)  are  the  even-  and  odd-  parts  of  the 
le      lo      2e      2o 

DHTs  H  (k)  and  H2(k). 

Expanding  the  right  hand  side,  we  get 

X*(k)X„(k)  =  [H,  (k)H.  (k)  +  H.  (k)H.  (k)] 
12        le     2e        lo     2o 

-  jfH,  (k)H.  (k)  -  B,  (k)H.  (k)]. 
le    2o       lo    2e 

Using  the  result  of  (3.16),  we  obtain  the  DHT  of  the  correlation  as 

Hfcross-correlation  of  x  (n)  and  x_(n)] 

=  H,  (k)H„  (k)  +  H,  (k)H.  (k)  +  H.  (k)H.  (k)  -  H,  (k)H.  (k) . 
le    2e    y     lo    2o       le    2o       lo    2e 

Substituting  for  H,  (k) ,  H,  (k),  H.  (k),  H.  (k)  in  terms  of  H  (k), 
le       lo      2e      2o  1 

D  (-k),  H  (k),  H  (-k)  and  then  simplifying,  we  get 

X  L,  /* 
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Hfcross  correlation  of  x  (n)  and  x.(n)] 

=  (l/2)[H1(k)E2(k)  +  H  (-k)H2(k)  -  H1(k)H2(-k)  +  H  (-k)H  (-k)]   (3.20) 

Symmetry  in  one  or  both  of  the  sequences  correlated  simplifies  the  above 
result.   Following  is  a  summary: 


SYMMETRY 


DHT  OF  THE  CROSS  CORRELATION 


{x  (n)}  is  even 


{x_(n)}  is  even 


Both  sequences  are  even 


{x  (n) }  is  odd 


{x  (n)}  is  odd 


Both  sequences  are  odd 


H 

(k)n2(k) 

H 

(-k)H2(k) 

H 

(k)n2(k) 

-H 

(k)H2(-k) 

H 

(k)n2(k) 

H 

(k)H2(k) 

Vfhen  we  need  to  perform  linear  cross-correlation  of  two  data  sequences, 
{x1(n)}  of  length  N  and  {x_(n)J  of  length  N  ,  we  add  a  string  of  zeros 

to  the  data  sequences  until  the  length  of  each  becomes  atleast 

(N  +  N.  -  1),  and  then  perform  a  circular  cross  correlation  as  shown 

above. 


3.4.11     Autocorrelation  Theorem 

Given    an    N-point    sequence     (x4(n)}     its    autocorrelation    R         is 

1  xx 

defined  as 

N-l 

Rxx(n)    =    (1/N>      ^     x    (h)x    (n  +   h).    n  =   0,    1.    ....    N-l. 
h=0 
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Observe  that  if  x  (n)  =  x_(n)  in  (3.19),  we  get  the  above  equation  for 

autocorrelation.   Therefore,  taking  the  result  for  the  DHT  of  the 
cross-correlation  for  two  sequences  in  (2.20),  and  making  the  substitu- 
tion H  (k)  ■  H  (k),  we  obtain  the  DHT  of  autocorrelation  of 

U^n)}. 

Therefore, 

R   (k)  =  (l/2)[H?(k)  +  H?(-k)].  (3.21) 

xx  11 

From  the  above  equation  it  is  clear  that  the  DHT  of  the  autocorrelation 

of  a  real  sequence  is  a  non-negative  even  function.   From  the  properties 

of  the  DHT  we  know  that  if  x  (n)  is  even,  H  (k)  =  E  (-k)  and  if  x  (n)  is 

7 
odd,  H  (k)  =  -  H  (-k).   In  either  case  (3.21)  reduces  to  H  (k) . 

If  we  need  to  perform  the  linear  autocorrelation  of  x  (n),  we  add  zeros 

to  the  data  sequence  until  its  length  becomes  atleast  2N  -  1  and  then 
perform  the  circular  autocorrelation  as  described  above  on  the  resulting 
sequence. 


3.4.12   Parseval 's  Theorem 

If  {x(p)}  is  an  N-point  real-valued  sequence,  we  know  that  its  DHT 

is  given  by 

N-l 
H(k)  =  (1/N)   ^  x(p)cas(2npk/N),  k  =  0,  1,  ...,  N-l.   (3.22) 
p=0 

By  changing  the  index  in  the  data  array  from  p  to  q,  we  can  again  write 

the  DHT  of  the  data  sequence  as 
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N-l 

H(k)  =  (1/N)   J   x(q)cas(2jTqk/N),  k  =  0,  1 N-  1.  (3.23) 

q=0 

Multiplying  (3.22)  and  (3.23),  we  get 

N-l   N-l 
H2(k)  =  (1/N)2  ^    ]>  x(p)x(q)cas(2npk/N)cas(2nqk/N). 
p=0   q=0 

k  =  0,  1.  ....  N  -  1. 

Now,  let  us  sun  both  sides  of  the  above  equation  over  N  points. 

N-l  N-l  N-l  N-l 

])  H^k)  =  (1/N)2   J   ^  x(P)x(^)   ^   cas(2npk/N)cas(2jiqk/N) 
k=0  p=0  q=0  k=0 

By  the  orthogonality  property  of  the  kernel,  we  have 

N-l 


y  cas(2npn/N)cas(2nqn/N)  =  N,  p  =  q 


n=0 

=  0,    otherwise 
Hence 

N-l  N-l 


]>     ^(k)    =    (1/N)        ]>     x2(p). 
k=0  p=0 


Changing  the  index  on  the  right  hand  side,  we  get 
N-l  N-l 


2  ^(k)  =  (1/N)   ^  x2(n). 


k=0  n=0 

The   above    equation   is   the  Parseval's    theorem    for    the    DHT.      Table    3.1 
provides   a    summary  of   the   properties   of   the   DHT. 
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H^k) 

H2(k) 

H^k)    - 

Hl(- 

-k) 

TABLE  3.1     PROPERTIES  OF  THE  DBT 

1.  Kernel   a   cas(2iran/N)  Orthogonal   and  periodic   on  N 

2.  (Effect   of  periodicity  (Imposes  periodicity   on    the 
of   the   kernel)                                                            data    sequence    and    its   DHT) 

3.  DHT{x(n)}  -  (Re{X(k)}   -   Im{X(k)}.    where 

X(k)    is   the  DFT  of  x(n)) 

4.  DFTU(n)}  -  (H    (k)   -  jH    (k),    where    H(k) 

e  o 

is   the  DHT  of  z(n)) 

5.  xx(n) 

6.  x2(n) 

7.  xx(n)    =  x1(-n) 

8.  xx(n)   -  -  xt(-n)  H^k)   -  -  H^-k) 

9.  xx(n)    +  x2(n)  H^k)    +  H2(k) 

10.  xx(n  +  m)  Ccos(2jnak/N)H1(k)   -   sin(2nmk/N)H1(-k)  ] 

11 .  [xx  (n)«x2 (n)  ]    (1/2)  [H1  (k)H2  (kJ+R^  <-k)H2  (k^  (k)H2 (-k)-^  (-k)H2 (-k)  ] 

12.  xx(n)x2(n)  (1/2)[H   (k)«H2(k)   +  H   <-k)«n2<k)    +  E^k^C-k) 

+  B   (-k)*H   (-k)] 

13.  {Cross-correlation  (1/2)  [H   (k)H   (k)   +  H1(-k)H2(k) 
of  xx(n)   and  x2(n))  -  H1(k)H2(-k)   +  B^-k^C-k)] 

14.  R        (x,(n)}  (l/2)[E*(k)   +  E*(-k)] 

xx       1  11 

N-l  N-l 

15.  (1/N)      ]>     x2(n)  -  ]>     B^k) 

n=0  k=0 
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4.   RADIX-2  DECIMATION- IN-TIME  FAST  HARTLEY  TRANSFORM  ALGORITHM 

4.1   Introduction 

We  know  that  the  DHT  of  an  N-point  real  sequence  {x(n)}  is  given  by 

N-l 
H(k)  =  (1/M)  ^  x(n)cas(2nnk/N),   k  =  0,  1,  ....  N-l. 
n=0 

However,     if  we    compute    the   DIIT  of    the    sequence    directly,    we    need    N 

multiplications    and    N  additions    to   compute    each   coefficient    in   the    DHT 

sequence,  and  a  total  of  w  multiplications  and  additions  to  compute  the 
N-point  sequence.  This  is  an  extremely  tall  order  for  long  data 
sequences.  As  such,  the  fast  algorithms  are  a  practical  way  out  of  this 
computational  explosion.  Bracewell  [7]  was  the  first  to  introduce  a 
radix-2  decimation- in-time  (DIT)  fast  algorithm  to  compute  the  DHT,  but 
he  did  not  exploit  many  inherent  symmetries  in  the  algorithm.  Kwong  and 
Shiu  [8]  proposed  a  more  refined  version  of  Bracewell 's  algorithm,  which 
significantly  reduced  the  number  of  multiplications  and  additions  re- 
quired to  compute  the  DHT.  Sorensen  e_t  a_l.  [9]  derived  the  same 
algorithm  by  an  index  mapping  approach.  This  algorithm  is  the  most 
commonly  used,  and  also  computationally  the  least  complex  of  all  the  FHT 
algorithms. 

However,  the  radix-2  algorithm  can  only  be  used  with  data  sequences 
whose  length  is  an  integer  power  of  two,  because  the  algorithm  involves 
successively  splitting  the  data  sequence  into  two  equal  half-length 
sequences.       Any    fast    algorithm   that    computes    the    DHT   of   a    given    length 
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of  data  sequence  from  smaller  length  DHTs  by  splitting  the  data  into 
smaller    sequences    is   called   a  DIT   algorithm. 

4.2     The   Decomposition  Formula 

By  exploiting  the  symmetry  and  the  periodicity  of  the  kernel 
cas(2nnk/N),  the  DHT  computation  can  be  decomposed  into  successively 
smaller  DHT  computations. 

Let    us    consider    the    DHT   equation    of    an    N-point    real-valued 

sequence. 

N-l 
H(k)    =    (1/N)      5     x<n)cas(2nnk/N).    k  =  0,    1,    ....    N-l.         (4.1) 
n=0 

Since  N  is  an  even  integer  (which  is  a  essential  for  this  algorithm),  we 
can  compute  H(k)  by  separating  {x(n)}  into  two  (N/2)-point  sequences, 
one  consisting  of  the  even-indexed  samples  and  the  other  containing  odd- 
indexed  samples  in  (x(n)}.   Stated  mathematically, 

t(k)  =   3   x(n)cas(2rtnk/N)  +  ^  x(n)cas(2nnk/N) 


HI 

n  even  n  odd 


By  substituting  the  variables  n  =  2r  for  even  n  and  n  =  2r  +  1  for  odd 

n,  we  get 

(N/2)-l 
H(k)  =   ^   x(2r)cas{2n(2r)k/N} 
r=0 

(N/2)-l 
+   ^     x(2r  +  l)cas{2n(2r  +  l)k/N}.  (4.2) 

r=0 
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However, 

cas{2*(2r  +   l)k/N}    =   cas{2n(2r)k/N  +  2nk/N} 

=   cos{2n(2r)k/N  +  2nk/N}    +   sin{2n(2r)k/N  +  2nk/N} 

=  cos{2n(2r)k/N}cos(2nk/N) 

-   sin{2n(2r)k/N}sin(2nk/N) 

+   sin{2ji(2r)k/N)cos(27ik/N) 

+   cos{27t(2r)k/N}sin(2:rk/N). 
Rearranging    the    terms,    we    get 

cas{2rt(2r  +  l)k/N}    =   cosUnk/N)    {cos[2n  (2r)k/N]    +   sin[2n(2r)k/N] } 

+   sin(2nk/N)    {cos [2n(2k)r/N]   -    sin[2n(2k)r/N] } . 

=   cos(27Tk/N)cas{2jr(2r)k/N} 

+   sin(2nk/N)cas{27r(2k)(-r)/N} 

=  cos(2nk/N)cas{2nrk/(N/2)} 

+  sin(2nk/N)cas{2nk(-r)/(N/2)} 

Substituting   the   above   expression   for   cas{2n(2r    +    l)k/N}    in    (4.2),    we 

get 

(N/2-1) 
H(k)  =   ]>  x(2r)cas{2nrk/(N/2)} 
r=0 

(N/2)-l 
+  cos(2nk/N)    J  *<2r  +  l)cas{2nrk/ (N/2) } . 
r=0 

(N/2)-l 
+  sin(2nk/N)    J  x(2r  +  Dca${2nk(-r)  /  (N/2) }  . 
r=0 

Recognizing  that  each  of  the  above  summations  is  an  (N/2)-point  DlfT,  we 

can  rewrite  the  above  equation  as 
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H(k)  =  H2r(k)  +  cos(27ik/N)H2r+1(k)  +  sin(2nk/N)H2r+1 (-k) 

=  H.  (k)  +  cos(27ik/N)H/.    (k)  +  sin(2nk/N)H„  ^{(N/2)  -  k} , 
£.1  zr+l  2r+l 

k  -  0,  1,  ....  (N  -  1).  (4.3) 

where 

H   (k)  is  the  DHT  of  the  (N/2)-point  even-indexed  samples  and 

£  T 

H„    _(k)    is   the  DHT  of   the    (N/2)-point    odd-indexed    samples.       Since    all 

Zr+l 

of   the   above   are    (N/2)-point   DHTs,    they   are   all   periodic   on  N/2. 

Therefore,    the    steps    involved    in  the    computation  of   an  N-point    DHT 
are : 

1.  Split  the  data  sequence  into  two  (N/2)-point  sequences,  one 
consisting  only  of  the  even-indexed  samples  and  the  other  consisting 
only  of   the   odd-indexed   samples. 

2.  Compute  the  DHTs  of  the  two  sequences  separately.  Call  the  DHT  of 
the    even-indexed    samples    H       (k)    and    the    DHT    of    the    odd-indexed 

samples   E.      ..  (k) . 

3.  Compute  the  sum  in  (4.3)  to  obtain  the  DHT  of  the  N-point 
sequence . 

However  we  can  repeatedly  apply  Step  1  in  performing  Step  2,  dividing 
each  (N/2)-point  sequence  into  two  (N/4)-point  sequences,  and  each 
(N/4)-point  sequence  into  two  (N/8)-point  sequences  and  so  on,  until  we 
reach  the  stage  where  we  need  to  compute  only  2-point  DHTs,  which  re- 
quire no  multiplications.  In  the  next  section,  the  working  of  the 
algorithm    is    explained  with  application   to    a   16-point   DHT. 
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4.3      Example   Description   of    the    Algoritlun 

Let  us  consider  a  16-point  data  sequence  {x(0),  x(l),  x(2),  .... 
x(15)}.  Now,  as  Step  1  suggests  in  Sec.  4.2,  let  us  rearrange  the  data 
points  such  that  the  first  half  of  the  array  consists  of  the 
even-indexed  samples  and  the  second  half  consists  of  the  odd-indexed 
samples.  Before  computing  the  DHTs  of  the  two  8-point  sequences  above, 
let  us  continue  splitting  them,  as  in  Step  1,  until  we  reach  the  2-point 
sequences.      Table   4.1    on   the    next    page    explains    the    procedure. 

The  last  column  in  Table  4.1  is  the  result  of  repeatedly  applying 
Step  1  to  the  data  sequence.  If  we  observe  closely,  the  index  of  any 
sample  in  the  last  column  is  obtained  by  taking  the  binary  repre- 
sentation of  the  data  sample  index  in  the  first  column  of  the  same  row, 
and  then  reversing  the  bit  string.  For  example,  let  us  consider  x(12) 
in  the  first  column.  The  number  12  is  represented  as  1100  in  binary 
form.  If  we  reverse  the  digits,  we  obtain  0011,  whose  decimal  equiv- 
alent is  3,  and,  as  such,  x(3)  is  the  sample  point  found  in  the  last 
column  for  the  row  containing  x(12).  Similarly,  the  number  14  is  repre- 
sented as  1110  in  binary  form,  and,  reversing  it  we  obtain  0111,  which 
is  the  binary  equivalent  of  the  number  7.  That  explains  why  we  find 
x(14)  in  the  first  column,  and  x(7)  in  the  last  column  of  the  same  row 
in  the   above    table. 

Therefore,  in  the  first  step  we  scramble  the  data  array  such  that 
it  contains  the  samples  in  the  bit  reversed  order.  The  scrambled  data 
is    shown   as    the    input    to    the    algorithm    in  Fig.    (4.1). 

From  Eqn.    (4.1),    we    see    that    the    DHT  of    a    two-point    sequence 
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(x(0),  x(l)}  is  obtained  by  once  adding  them,  and  next  subtracting  one 
from  the  other    i.e., 

H(0)  =  x(0)  +  x(l). 
and 

H(l)  =  x(0)  -  x(l). 

TABLE  4.1   SCRAMBLING  THE  DATA  FOR  THE  RADIX-2  DIT  ALGORITHM 


Original 

1st 

2nd 

3rd 

array 

scramble 

scramble 

scramble 

x(0) 

x(0) 

x(0) 

x(0) 

x(l) 

x(2) 

x(4) 

x(8) 

x(2) 

x(4) 

x(8) 

x(4) 

x(3) 

x(6) 

x(12) 

x(12) 

x(4) 

x(8) 

x(2) 

x(2) 

x(5) 

x(10) 

x(6) 

x<10) 

x(6) 

x(12) 

x(10) 

x(6) 

x(7) 

x(14) 

x(14) 

x(14) 

x(8) 

x(l) 

x(l) 

x(l) 

x(9) 

x(3) 

x(5) 

x<9) 

x(10) 

x(5) 

x(9) 

x(5) 

xUl) 

x(7) 

x(13) 

x<13) 

x(12) 

x(9) 

x(3) 

x(3) 

x(13) 

x(ll) 

x(7) 

x(H) 

x(14) 

x(13) 

x(ll) 

x(7) 

x(15) 

x(15) 

x(15) 

x(15) 
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Computation  of  the  2-point  DHTs  is  shown  in  Fig.  (4.1).  The  broken 
lines  in  the  figure  indicate  subtraction  and  the  solid  lines  indicate 
addition.  From  these  two-point  DIITs,  we  compute  the  four-point  DHTs  in 
the  second  stage.  Computation  of  a  four-point  DIIT  also  does  not  need 
any  multiplications. 

In  the  third  stage,  we  compute  the  8-point  DIITs  from  the  4-point 
DHTs  using  (4.3),  bearing  in  mind  that  the  top  4-point  sequence  always 
acts   as  H      (k)    and   the   bottom   one    as   H  (k).      The    internal    structure 

of    the    boxes  marked  T      in  Fig.    (4.1)    is    shown   separately    in  Fig.    (4.2). 

The    input    to    the    bozes    labeled   T       is    always    the    bottom    sequence 

H  (k).       The    boxes    implement    the    summation    cos(2nk/N)H  (k)    + 

sin(2nk/N)H„         f(N/2)-k)   with  N  =   8,    k  =   0.    1,    ....    3. 

Zr+1 

Finally  the  16-point  sequence  is  obtained  from  the  two  8-point 
sequences,    with   the    top   8-point    sequence    acting   as   H      (k)    and    the   bottom 

8-point    sequence    acting    as    H„       ,  (k)    in  Eqn.    (4.3).      The   box  T.    imple- 

2r+l  4 

ments  the  summation 

cos(2jik/N)H<.  -  (k)  +  sin(2nk/N)H.   ,{(N/2)-k},  N  =  16,  k  =  0,  ....  7 
2r+l  2r+l 


4.4   Conmvtational  Cost 

Let  us  consider  Eqn.  (4.3). 

H(k)  =  n.  (k)  +  cos(2nk/N)n.  _ (k)  +  sin(2nk/N)B.  ^,{(N/2)  -  k} 
lt  zr+i  zr+l 

From  the  above  equation  we  can  form  the  following  three  equations. 
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HUN/2)     -    k}     =    H.     {(N/2)     -    k}     -    cos(2jtk/N)H,        ,{(N/2)     -     k} 

It  2  r  +  1 


+   sin(2nk/N)H2r+1(k)  (4.4) 


H(k  +  N/2)    =  H„    (k) 
2r 


-    [cos(2nk/N)H_         (k)    +   sin(2nk/N)H.    ^{(N/2)   -  k}],      (4.5) 
Zr+1  2r+l 

H(N  -   k)        =  H„    {(N/2)    -   k} 
2r 

+    [cos(2nk/N)H.    ^{(N/2)   -  k}    -   sin(2nk/N)H„      ,  <k) .         (4.6) 
2r+l  2r+l 

In  arriving   at   the   above   equations,    we    used    the    fact    that    H      (k)    and 

H       ..  (k)    are   periodic   on  N/2. 

By  referring   to    (4.3),    (4.4),    (4.5),     (4.6)    we    see    that    the    basic 

products    involved    in    those    equations    are    only   those    that    are    shown   in 

Table  4.2 

TAELE  4.2   PRODUCTS  REQUIRED  TO  COMPUTE  AN  N-POINT  DHT 

USING  THE  RADIX-2  DIT  ALGORITHM 

PI  =  cos(2nk/N)[H-  ^(k)]    P2  =  cos  (2nk/N)  [H.  ^{(N/2)  -  k}  ] 
zr+i  zr+i 

P3  =  sin(2nk/N)[H.  ^  (k)]    P4  =  sin(2nk/N) [H,  ^{(N/2)  -  k} ] 
2r+l  2r+l 

By  forming  the  above  products,  we  observe  that  for  each  k  we  can  compute 
a  total  of  four  points  n(k),  H{(N/2)  -  k}  .  H{(N/2)  +  k}  and  H(N  -  k) . 
Thus,  once  we  form  the  products  necessary  to  compute  H(0),  H(l),  .... 
H{(N/4)},  we  can  compute  the  remaining  (3N/4)  -  1  coefficients  in  the 
DHT  sequence  without  any  further  multiplications.  Table  4.3  explains 
the  fact  clearly  for  a  16-point  sequence. 
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TABLE  4.3      WORKING  OF  THE   RAIUX-2   DIT  AGORITHM      FOP.  A   16-POINT  DHT 

N  -   k 


15 
14 
13 


k 

.(N/2)    - 

k 

(N/2)    +  k 

0 

8 

1 

7 

9 

2 

6 

10 

3 

5 

11 

4 

12 

Thus,  for  a  16-point  data  we  can  compute  the  entire  DHT  sequence  by 
forming  the  following  products 

PI,  P2,  P3,  P4,    k  =  0,  1,  ....  4 
in  Table  4.2,  and  then  using  them  in  Eqns.  (4.4)  -  (4.6),  instead  of 
directly  computing  the  coefficients  using  Eqn.  (4.3). 

From  Table  4.2  it  is  clear  that  whenever  k  =  0  or  k  =  N/4,  the 
trigonometric  terms  in  products  PI  through  P4  become  either  zero  or  one 
and  as  such,  we  do  not  need  any  multiplications  for  those  two  special 
cases.  Avoiding  these  two  special  cases,  we  need  to  compute  the  four 
products 

PI,  P2,  P3.  P4.    k  =  0,  1,  ...,  (N/4)  -  1. 
So  for  an  N-point  sequence  we  need  4[(N/4)  -  1]  =  N  -  4  products. 

Since  we  repeatedly  apply  Eqn.  (4.3)  in  the  algorithm  by  computing 
(N/2)-point  DHTs  from  two  (N/4)-point  DHTs,  and  (N/4)-point  DHTs  from 
two  (N/8)-point  DITTs  etc.,  we  can  reap  the  above  mentioned  advantages  in 
multiplications  at  every  stage.  Therefore  the  total  number  of  multi- 
plications is  given  by 
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M  =  (N  -  4)  +  2[(N/2)  -  4]  +  22[(N/22)  -  4]  +  ...  +  2?  3[N/(2?  3)  -  4] 

where 

P  =  log2N. 

Simplifying  the  above  expression,  we  get 

M  =  (P  -  2)N  -  4(1  +  2  +  22  +  23+  ...  +  2P~3) 

=  (P  -  2)N  -  4(2P"2  -  1) 
=  (P  -  2)N  -  N  +  4 
=  PN  -  3N  +  4. 
Hence  total  number  of  multiplications  is  given  by 

M  =  N  log  N  -  3N  +  4. 

Additions 

Having  formed  the  products 

PI.  P2,  P3,  P4,    k  =  0.  1,  ....  (N/4)  -  1, 
let  us  form  the  following  two  summations 

51  =  PI  +  P4   and 

52  =  P3  -  P2 

Then   (4.3),    (4.4),    (4.5),    (4.6)    can  be    rewritten  as 
H(k)   =  H      <k>    +   S1, 
H{(N/2)    +  k}    =   H2j.(k)    -    SI, 

H{(N/2)    -  k}    =  H„    {(N/2)    -   k}    -   S2, 
2r 

H(N  -  k)      =  H„    {(N/2)   -  k}    +   S2. 
Zr 
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Therefore,  in  addition  to  the  two  summations  SI  and  S2  ,  we  need  the 

above  four  summations,  with  k  taking  values  from  1  through  (N/4)  -  1), 

to  compute  the  entire  N-point  DKT  sequence.   Thus  a  total  of 

6[(N/4)  -  1]  =  (3N/2)  -  6  additions  are  needed.   But  for  the  two  special 

cases  of  k  =  0  and  k  =  N/4,  we  need  only  two  additions  each  instead  of 

the  usual  six   as  explained  below. 

Substituting  k  =  0  and  k  =  (N/4)  in  (4.3),  (4.4),  (4.5)  and  (4.6),  we 

get 

11(0)  =  B.  (0)  +  H0  ^(0) 
lt  zr+l 


and 


H(N/2)  =  B„  (0)  -  B.  ^(0). 
2r       2r+l 

H(N/4)  =  B.  (N/4)  +  B.  ^(N/4), 

and 

B(3N/4)  =  B.  (N/4)  -  B.  _ (N/4) . 
zr         zr+i 

Therefore,  total  number  of  additions  A  is  given  by 

A  =  (3N/2)  -  6  +  4  =  (3N/2)  -  2. 

Since   we   continue    splitting   the   data    sequence    into  half-length   sequences 

until    we    reach   2-point    sequences,    we    can   obtain   the   above    savings    in 

additions    at    every    stage.      Therefore,    total    number    of    additions    A 

is   given  by 

A  =  [(3N/2)  -  2]  +  2[(3/2)(N/2)  -  2]  +  ...  +  2?~2[ (3/2) (N/2P"2)  -  2] 
where 

P  =  log2N. 

A  =  (3N/2HP  -  1)  -  2[1  +  2  +  22+  ...  +  2P~2] 
=  (3N/2MP  -  1)  -  2[2P~1  -  1] 


=    (3N/2XP  -   1)   -   N  +  2. 
The   above    is    the    number    of    additions    needed    until    we    reach   4-point 
sequences.      In  addition,    we   have    (N/2)    two-point    sequences,    which   need  N 
more   additions 

A  =    (3N/2)P  -    (3N/2)    -   N  +  2   +  N 

=    (3N/2)log2N  -    (3N/2)    +  2. 
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5.   RADIX-2  DECIMATION- IN-FREQDENCY  FAST  HARTLEY  TRANSFORM  ALGORITHM 

5.1  Introduction 

Decimation-in-time  (DIT)  and  decimat  ion- in-f  requency  (DIF)  algo- 
rithms differ  in  one  significant  aspect.  DIT  algorithms  compute  the 
DHTs  from  short  length  DIITs  and  DIF  algorithms  compute  the  DIITs  a_s  short 
length  DHTs.  That  is,  when  we  compute  the  DHT  of  a  16-point  sequence 
using  DIT  algorithm,  we  first  compute  the  2-point,  4-point,  8-point  DHTs 
and  finally  by  combining  the  two  8-point  DIITs,  we  obtain  the  16-point 
DHT,  In  a  DIF  algorithm,  the  task  of  computing  a  16-point  DHT  is 
progrssively  reduced  to  that  of  computing  only  2-point  DHTs,  without 
ever  actually  computing  the  intermediate  8-point  and  4-point  DHTs. 
Computationally,  a  DIF  algorithm  does  not  have  any  advantage  over  the 
corresponding  DIT  algorithm.  A  radix-2  DIF  algorithm  can  be  used  only 
with  data  sequences  whose  length  is  an  integral  power  of  two. 
Meckelburg  and  Lipka  [10]  first  introduced  the  radix-2  DIF  FHT  algo- 
rithm, Sorensen  e_t  a_l.  [9]  derived  the  same  algorithm  by  an  index 
mapping   approach. 

5.2  The  Decomposition  Formula 

The   DHT  H(k)    of    an  N-point    real    data    sequence    {x(n)}    is    given  by 

N-l 
H(k)    =    (1/N)      J     x(n)cas(2nnk/N),    k  =   0.    1,    ...,    N-l.  (5.1) 

n=0 

We    can  divide    {x(n)}    into    two   half-length    sequences    and    rearrange 

the  above  equation  as 
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(N/2)-l 
H(k)  =  (1/N)    ^    x(n)cas(2nnk/N) 
n=0 

(N/2)-l 
+  (1/N)    ^    x(n)cas(2nnk/N).  k  =  0,  1.  ...,  N  -  1. 
n=N/2 


or 

(N/2)-l 
H(k)  =  (1/N)    ^    x(n)cas(2nnk/N) 
n=0 

(N/2)-l 
+  (1/N)    ^    x(n  +  N/2)cas{2nk(n  +  N/2)/N),  k  =  0,  1,  ...,  N  -  1. 
n=0 

(N/2)-l 
H(k)  =  (1/N)    }   x(n)cas(2nnk/N) 
n=0 

(N/2)-l 
+  (1/N)    ^   x(n  +  N/2)cas[kn  +  2nnk/N) ,  k  =  0,  1,  ...,  N  -  1.(5.2) 
n=0 

Now,  let  us  compute  the  even-  and  odd-  indexed  DHT  samples  separately. 

Substituting  k  =  2r  in  (5.2),  we  obtain  the  even-indexed  samples 

(N/2)-l 
H(2r)  =  (1/N)    ^   x(n)cas{2nn(2r)/N} 
n=0 

(N/2)-l 
+  (1/N)    ^  x(n  +  N/2)cas[2rjr  +  2nn(2r)/N],  r  =  0.  1,  ...,  (N/2)  -  1. 
n=0 

Recognizing  that 

cas[2rn  +  2nn(2r)/N]  =  cas(2nn(2r)/N) 

we  can  rewrite  the  above  equation  as 
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(N/2)-l 
H(2r)  -  (1/N)    J    [x(n)  +  x(n  +  N/2)]cas{2nn(2r)/N)} 
n=0 

(N/2)-l 
=  (1/N)    ^   [x(n)  +  x(n  +  N/2) ]cas [2nnr/ (N/2) ] , 
n=0 

r  =  0,  1.  ....  (N/2)  -  1.       (5.3) 

Since  the  above  is  the  expression  for  an  (N/2)-point  DHT,  we  conclude 

that  the  even  indexed  samples  of  the  DHT  can  be  computed  as  an 

(N/2)-point  DHT. 

Substituting  k  =  2r  +  1  in  (5.2),  we  obtain  the  odd-indexed  samples 

of  the  DHT  sequence. 

(N/2)-l 
H(2r  +  1)  =  (1/N)    ^   x(n)cas[2nn(2r  +  1)/N] 

n=0 

(N/2)-l 
+    (1/N)        ^        x(n  +  N/2)cas[(2r  +   l)n   +  2nn(2r  +   1)/N].       (5.4) 
n=0 

Using   the    fact   that 

cas[(2r  +   l)n   +  2nn(2r  +   1)/N]    =   cas[n   +  2nn(2r  +   1)/N] 

=  -   cas[2nn(2r  +   1)/N] 

we   can  rewrite    (5.4)   as 

(N/2)-l 
H(2r  +   1)    =    (1/N)        ^       x(n)cas[2nn(2r  +   D/N] 

n=0 

(N/2)-l 
-    (1/N)        ^       x(n  +  N/2)cas[2nn(2r  +   D/N], 
n=0 

r  =   0,    1.     ....    (N/2)    -   1.         (5.5) 
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(N/2)-l 

=    (1/N)        2        tx(n)   "  x(n  +  N/2)]cas[2nn(2r  +   1)/N].  (5.6) 

n=0 

However 

cas[2rrn(2r  +  1)/N]    =  cos (2nn/N)cas(2nn(2r)/N) 

+   sin(2nn/N)cas[2n(2r)(-n)/N] 

Substituting   the   above   result    in   (5.6),    we   get 

(N/2)-l 
H(2r   +   1)    =    (1/N)        ^        [x(n)    "   x(n  +  N/2)]cos(27tn/N)cas(2Ttnr/(N/2) 

n=0 

(N/2)-l 
+    (1/N)      )        [x(n)   -  x(n+N/2)]sin(2nn/N)cas[2nr(-n)/(N/2)]    (5.7) 
n=0 

Reversing  the  direction  of  the  second  summation  in  (5.7),  we  get 

(N/2)-  1 
H(2r+1)  =  (1/N)    5   C[x(n)  ~  *<n  +  N/2)]cos(2nn/N) 
n=0 

+  [x{(N/2)  -  n}  -  x(N  -  n)]sin(2Trn/N)}cas{2nnr/(N/2)}     (5.8) 

The  above  is  an  (N/2)-point  DHT  to  compute  the  odd-indexed  samples  of 

the  DHT. 

Hence  from  (5.3)  and  (5.8)  we  see  that  an  N-point  DHT  can  be  computed  as 

two  (N/2)-point  DHTs.   The  following  are  the  steps  involved  in  the 

computation  of  the  DHT  of  an  N-point  sequence  (x(n)}. 

1.  First,  form  the  N/2-point  sequences  {g(n)}  and  (h(n)},  where  g(n) 
and  h(n)  are  given  by 

g(n)  =  [x(n)  +  x(n  +  N/2)],  n  =  0,  1,  ....  (N/2)  -  1, 
h(n)  =  [x(n)  -  x(n  +  N/2)],  n  =  0,  1.  ....  (N/2)  -  1. 

2.  Then  compute  the  (N/2)-point  sequence 
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f(n)    =   h(n)cos(2nn/N)    +   h{(N/2)    -   n} sin(2nn/N) , 

n  =   0,    1 (N/2)    -    1. 

3.      Compute    the    (N/2)-point    DHTs   of    the    sequences    g(n)    and    f(n). 

The   DHT  of   g(n)    gives    the    even- indexed    samples,    and    the   DHT    of    f(n) 

gives  the  odd-indexed  samples  of  Il(k). 
However,  in  performing  Step  3  we  can  apply  Steps  1  and  2  to  the  computa- 
tion of  those  (N/2)-point  DHTs  repeatedly,  until  we  need  to  compute  only 
2-point  DHTs.  It  is  this  repeated  application  of  Steps  1  and  2  that 
makes  the  computation  of  the  DHT  faster.  The  working  of  the  algorithm 
for   an  8-point   data    sequence    is   explained   in  the   next    section. 

5.3      Example   Description  of   the   Algorithm 

Fig.  (5.1)  shows  the  flow  diagram  for  an  8-point  FHT.  Observe  that 
in  Fig.  (5.1),  the  data  on  the  input  side  of  the  flow  diagram  is  in 
the  normal  order.  We  do  not  need  to  scramble  the  data  in  the  beginning 
as  we  did  with  the  DIT  radix-2  algorithm.  In  the  diagram,  whenevr  two 
solid  lines  meet  at  a  point,  it  corresponds  to  an  addition  and  when  a 
solid   line   and   a  broken   line  meet,    it   defines    subtraction. 

As  suggested  in  Step  1  in  the  previous  section,  we  first  form  the 
4-point    sequences    {g(n)}    and    {h(n)}    as    follows: 

g(n)    =   x(n)    +  x(n  +  4),    n  =   0,    1,    ...,    3 

h(n)  =  x(n)  -  x(n  +4),  n  =  0,  1,  ....  3 
The  box  marked  T  achieves  the  computation  of  f(n).   The  internal  struc- 
ture of  T.  is  shown  separately.   The  box  T   takes  the  sequence  h(n)  as 
its  input  and  computes  f(n)  as 
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f(n)  =  h(n)cos(2nn/N)  +  h(4  -  n)sin(2nn/N) ,  N  =  8,  n  =  0,  ...»  3. 
According  to  Step  3  in  the  previous  section,  we  now  compute  the 
DHTs  of  g(n)  and  f(n).  But  we  can  apply  Steps  1  and  2  to  those  two 
sequences  one  more  time,  leading  us  to  the  computation  of  only  2-point 
DHTs.  We  observe  that  the  DIIT  coefficients  are  in  the  bit  reversed 
order.  We  have  to  perform  a  bit  reversal  operation,  as  described  in  the 
previous  chapter,  to  get  the  coefficients  back  in  the  normal  order.  In 
this  example,  we  have  only  three  stages  involved,  reduction  of  an 
8-point  computation  to  two  4-point  DHTs,  and  4-point  DHTs  to  2-point 
DHTs   and   finally   the   computation   of    2-point    DHTs.       In    general,    for    a 

p 
data    sequence    of    length   N   =    2    ,    we    will    have   P   iterations    in  the  DHT 

computation  process.  Like  the  DIT  radix-2  algorithm  in  the  last  chap- 
ter, this  DIF  algorithm  also  is  an  in-place  algorithm.  Therefore,  once 
we  compute  g(n)  and  f(n)  from  the  original  8-point  data  in  the  first 
iteration,  we  can  replace  the  data  elements  with  g(n)  and  f(n)  in  the 
data  array.  We  keep  doing  this  in  all  the  iterations  until  we  compute 
the  required  2-point  DHTs  in  the  final  iteration  at  which  point  the 
original  data  elements  in  the  array  will  have  been  replaced  by  the  DHT 
coefficients. 

5.4      Computational   Cost 

We    stated   earlier   that    from   the  N  data   elements    in    (x(n)},    we    form 
the    (N/2)-point    sequences    {g(n)}    and    {h(n)},    and   from  h(n)    the    sequence 
f(n).      The    sequence    f(n)    is   given  by 
f(n)    =   h(n)cos(2rrn/N)    +  h{(N/2)   -  n}sin(2nn/N) ,    n  =  0,    ....    (N/2)    -   1. 
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Thus,  to  obtain  a  point  f(n)  we  need  two  multiplications 

1.  h(n)cos(2nn/N) 

2.  h{(N/2)-n}sin(2nn/N). 

Since  {f(n)}  contains  a  total  of  (N/2)  points,  we  thus  need  a  total  of 
2(N/2)  =  N  multiplications.  However,  whenever  n  =  0  or  n  =  N/4,  the 
trigonometric  functions  in  the  above  two  multiplications  reduce  to 
either  zero  or  one,  and,  as  such,  we  can  avoid  those  multiplications. 
Therefore,  those  two  special  coefficients  save  us  a  total  of  four 
multiplications.  Hence,  the  number  of  multiplications  required  to 
obtain  the  (N/2)-point  sequence  {f(n)}  from  the  N-point  data  is  (N-4). 
These  savings  in  multiplications  can  be  obtained  in  all  the  iterations. 
Adding  all  such  multiplications  until  we  reach  4-point  DFITs  (wherefrom 
we  do  not  need  any  multiplications),  the  total  count  M  is  given  by 

M  =  (N-4)  +  2[(N/2)  -  4]  +  22[(N/22)  -  4]  +  ...  +  2P"3 [ (N/2P~3 )  -  4] 
where  P  =  log.N. 

Simplifying  further,  we  get 

M  =  (P  -  2)N  -  4[1  +  2  +  22  +  ...  +  2P~3] 

-  (P  -  2)N  -  4[2P"2-1] 

=  (P  -  2)N  +  4  -  N        (since  2P  =  N) 

=  PN  -  3N  +  4 
M  =  Nlog  N  -  3N  +  4. 

Let  us  consider  the  number  of  additions.   From  the  original  data  se- 
quence {x(n)},  we  first  form  g(n)  and  h(n)  as  shown  below. 
g(n)  =  tx(n)  +  x(N  +  N/2)],  n  =  0,  1.  ....  (N/2)  -  1 
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h(n)  =  [x(n)  -  x(n  +  N/2)],  n  =  0.  1.  ....  (N/2)  -  1. 

We  need  (N/2)  additions  to  form  {g(n)J  and  (N/2)  additions  to  form 
(h(n)},  totaling  N  additions.  Besides,  in  the  computation  of  f(n),  we 
come  across  additions  of  the  following  type: 

f(n)  =  h(n)cos(2nn/N)  +  h[(N/2)  -  n] sin(2rm/N) ,  n  =  0,  ....  (N/2)  -  1. 
We  thus  need  (N/2)  more  additions  to  compute  {f(n)}.  However,  when 
n  =  0  or  n  =  (N/4),  the  trigonometric  functions  in  the  above  addition 
become  either  zero  or  one  and,  as  such,  we  do  not  need  any  additions  in 
evaluating  f(0)  and  f(N/4),  thus  saving  us  two  additions  in  computing 
(f(n)}.  Therefore,  to  compute  {f(n)l  we  need  [(N/2)  -  2]  additions. 
Hence,  the  total  number  of  additions  in  computing  (g(n)},  (h(n)},  (f(n)} 
are  given  by 

Al  =  (N/2)  +  (N/2)  +  (N/2)  -  2 
=  (3N/2)  -  2 
Summing  such  additions  over  all  the  iterations  until  we  reach  2-point 
DHTs,  the  total  number  of  additions  A  is  given  by 
A  =  [(3N/2)  -  2]  +  2[(3/2)(N/2)  -  2] 

+  22[(3/2)(N/22)  -  2]  +  ...  +  2P"2[(3/2)(N/2P"2)  -  2]. 

=  (3N/2MP  -  1)  -  2[1  +  2  +  22  +  ...  +  2P"2] 

=  (3N/2MP  -  1)  -  2[2P~1  -  1] 

=  (3N/2MP  -  1)  -  N  +  2.       (since  2P=  N) 
In  addition,  we  will  have  (N/2)  two-point  DHTs  to  compute,  thus  requir- 
ing 2(N/2)  =  N  more  additions. 
Thus,  number  of  additions  is  given  by 
A  =  (3N/2MP  -D-N  +  2  +  N 
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-    (3N/2)P  -   3N/2    +   2 


=    (3N/2)log2N  -   3N/2    +  2. 


The    above    count    for    additions   and  multiplications    is    the    same    as    that 
for   the   DIT  radix-2   algorithm  discussed    in  the   previous   chapter. 
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6.   RADIX-4  DECIMATION- IN-TIME  FAST  HARTLEY  TRANSFORM  ALGORITHM 

6.1  Introduction 

In  the  previous  two  chapters,  we  discussed  the  radix-2  algorithms 
which  can  be  used  to  compute  the  DHT  of  a  sequence  whose  length  is  an 
integer  power  of  two.  In  addition  to  being  even,  if  the  length  of  a 
sequence  is  an  integer  power  of  four,  we  can  use  a  radix-4  algorithm  to 
compute  its  DHT  more  efficiently.  In  a  DIT  radix-4  algorithm,  we  split 
the  data  sequence  into  four  equal-length  smaller  data  sequences,  compute 
the  DIITs  of  those  sequences,  and  finally  combine  them  to  obtain  the  DHT 
of  the  original  data.  Sorensen  e_t  a_l.  [9]  derived  a  radix-4  algorithm 
by  an  index  mapping  approach.  Prado  [11]  also  presented  a  radix-4 
algorithm,  which  takes  more  number  of  operations  than  that  of  Sorensen 
et  a  1 .  In  this  chapter,  we  derive  the  radix-4  algorithm  from  the 
definition  of  the  DHT,  without  relying  on  the  index  mapping  approach, 
describe  its  working  with  an  example,  and  finally  obtain  its  computa- 
tional   cost. 

6.2  The  Decomposition  Formula 

The  DHT  of  an  N-point  real  data  sequence  (x(n)}  is  given  by 

N-l 
H(k)  =  (1/N)   5  x(n)cas(2nnk/N),   k  =  0.  1,  ....  N-l.  (6.1) 

n=0 

Since  the  length  of  (x(n)}  is  an  integer  power  of  four  (which  is  a 

precondition  for  using  this  algorithm),  we  can  split  the  summation  in 
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the  above  equation  into  four  equal-length  smaller  summations  as  shown 

below.   Since  the  factor  (1/N)  is  only  a  constant  outside  the  summation, 

we  need  not  carry  it  in  the  derivation  process  until  we  reach  the  final 

step. 

(N/4)-l  (N/4)-l 

H(k)    =       2  x(4n)cas[2n(4n)k/N]    +       J  x(4n  +   Dcas[2n(4n  +   l)k/N] 

n=0  n=0 

(N/4)-l 
+        5  x(4n  +  2)cas[2n(4n  +  2)k/N] 

n=0 

(N/4)-l 
+       ^  x(4n  +  3)cas[2n(4n  +  3)k/N].  (6.2) 

n=0 

Observing  that 

cas[2n(4n  +  l)k/N]  =  cos[2n(4n  +  l)k/N]  +  sin[2n(4n  +  l)k/N], 
and  then  using  the  familiar  trigonometric  identities,  we  get 

cas[2n(4n  +  l)k/N]  =  cos(2nlk/N)cas[2nnk/ (N/4) ] 

+  sin(2nlk/N)cas[27rn(-k)/(N/4)].  (6.3) 

Using  the  result  of  (6.3)  in  (6.2)  with  1  =  1,  2,  3,  we  get 

(N/4)-l 
H(k)  =   ^    x(4n)cas[2nnk/(N/4)] 
n=0 

(N/4)-l 
+  cos(2nk/N)    ^   x(4n  +  l)cas[2nnk/(N/4)] 
n=0 

(N/4)-l 
+  sin(2nk/N)   ^    x(4n  +  l)cas[2nn(-k)/(N/4) ] 
n=0 
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(N/4)-l 
+  cos[2n(2k)/N]    J   x(4n  +  2)cas[2nnk/ (N/4) ] 
n=0 

(N/4)-l 
+  sin[2n(2k)/N]   ^   x(4n  +  2)cas [2nn(-k)/ (N/4) ] 
n=0 

(N/4)-l 
+  cos[2n(3k)/N]   ^   x(4n  +  3)cas[2nnk/(N/4) ] 
n=0 

(N/4)-l 
+  sin[2n(3k)/N]   ^   x(4n  +  3)cas[2nn(-k)/ (N/4) ] , 
n=0 

k  =  0,  1,  ....  (N/4)  -  1.    (6.4) 

Bringing  back  the  (1/N)  factor  and  recognizing  that 

(N/4)-l 
(1/N)    ^   x(4n  +  l)cas[2nn(-k)/(N/4)],   k  =  0,  1.  ...,  (N/4)  -  1 
n=0 

is   an    (N/4)-point   DHT,    and    indicating    it    as   H.      ,,    we    can   rewrite     (6.4) 

4n+l 

as 

H(k)    =    [H.    (k)    +  cos(2nk/N)H.       . (k)    +   sin(2nk/N)ir       . (-k) 
4n  4n+l  4n+l 

+   cos[2n(2k)/N)H.      .(k)    +    sin[2n(2k)/N)H.      _ (-k) 
4n+2  4n+2 

+  cos[27t(3k)/N)IT   ,(k)  +  sin[2n(3k)/N)H\   ,(-k)]   (6.5) 
4n+3  4n+3 

Eqn.  (6.5)  is  the  decomposition  formula  for  the  DIT  radix-4  algorithm. 

Thus,  H(k)  is  computed  from  four  (N/4)-point  DHTs  -   H,  (k).  H,   ,  (k), 

4n  4n+l 

H.     ,„(k)    and    H.     ,_(k).       The    following    are    the    steps    involved    in   the 
4n+2  4n+3 

algorithm. 
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1.  Split  the  data  sequence  {x(n)}  into  four  (N/4)-point  sequences 
{x(4n)},  {x(4n  +  1)},  {x(4n  +  2)},  {x(4n  +  3)}. 

n  =  0,  1,  ...,  (N/4)  -  1. 

2.  Compute  independently  the  (N/4)-point  DHTs 

*a    (k)'  Ba    xi(k)»  Ea   o.o(k)  and  Ea   ..<»<*>'  k  =  0,  1.  ...,  (N/4)  -  1 
4n      4n+l      4n+2         4n+3 

respectively,  of  the  sequences  formed  in  Step  1. 

3.  Using  (6.5)  combine  the  above  (N/4)-point  DHTs  to  obtain  the  N-point 
DHT. 

6.3   Example  Description  of  the  Working  of  the  Algorithm 

Let  us  consider  a  64-point  sequence  (x(n)}.   As  Step  1  suggests  in 
Sec.  6.2,  let  us  divide  it  into  four  16-point  sequences  as  shown  below. 

(x(4n)}    =  {x(0),  x(4),  x(8),  ...,  x(60)} 

{x(4n  +  1)}  =  {x(D.  x(5).  x(9).  ....  x(61)} 

{x(4n  +  2)}  =  {x(2),  x(6),  x(10) x(62)} 

{x(4n  +  3)}    =    {x(3),    x(7).    x(ll),    ....    x(63)}. 
Now,    according   to   Step  2,    we   have    to   compute    the   DHTs   of   the   above 
16-point    sequences.       But    in   doing   so,    we   can  again  apply  Step   1,    thus 
splitting    each    16-point    sequence    into    four   4-point    sequences.      The 
division  of    {x(4n)}    results    in  the   following: 

(x(16n)}  =  {x(0).  x(16),  x(32),  x(48)} 
{x(16n  +  4)}  =  U(4),  x(20),  x(36),  x(52)} 
{x(16n  +  8)}  =  {x(8),  x(24),  x(40),  x(56)} 
{x(16n  +   12)}    =      {x(12).    x(28),    x(44).    x(60)} 


Similarly  when  we    split    the    second    16-point    sequence     {x(4n    +    1)},    we 

obtain    the    following    four-point    sequences. 

{x(16n  +   1)}      =      {x(l).    x(17),    x(33),    x(49)} 

{x(16n   +   5)}      =      {x(5),    x(21).    x(37).    x(53)} 

{x(16n   +   9)}      =      {x(9),    x(25),    x(41).    x(57)} 

{x(16n  +   13)}    =      {x(13),x(29).    x(45),    x(61)}. 

When  v/e    split    the    third    16-point    sequence     {x(4n    +    2)}     into    4-point 

sequences,    we   obtain 

{x(16n   +  2)}      =         {x(2),    x(18),    x(34),    x(50)} 

{x(16n  +   6)}      =         {x(6),    x(22),    x(38),    x(54)} 

{x(16n  +   10)}    =        {x(10),    x(26),    x(42),    x(58)} 

{x(16n  +  14)}    =        (x(14),    x(30),    x(46),    x(62)} 

And    finally   the    division   of    the    last    16-point    sequence    {x(4n+3)}    results 

in  the    following   four    sequences. 

{x(16n  +  3)}      =  {x(3),    x(19)»    x(35),    x(51)} 

{x(16n  +  7)}      =  {x(7).    x(23),    x(39).    x(55)} 

{x(16n  +   11)}    =  {x(ll),    x(27).    x(43).    x(59)} 

{x(16n  +  15)}    =  {x(15),    x(31),    x(47),    x(63)}. 

Therefore,    the    first    thing   that   needs    to  be    done    is    to    scramble    the    data 

array    such    that    it    contains    the    data    elements    in    the    above    order   of 

rows.      Some   kind   of    scrambling   of    the    input    data    is    the    characteristic 

of   any  DIT   algorithm. 

Once  we  scrambled  the  data  array,  the  first  stage  in  the  computa- 
tion process  is  to  compute  the  DHTs  of  all  the  above  4-point  sequences, 
and   this    does   not    require    any  multiplications. 
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Now    that  we   have    all    the   necessary   sixteen  4-point   DITTs,    we    coinbine 

2 
them  using    (6.5)    with  N  =   4      and   k  =   0,     1,     ...,     15    to    obtain    the    four 

16-point   DHTs    in   the    second    stage. 

Once   we   have    the    four   16-point   DHTs,    we    compute    the    64-point    DHT 

3 
from    them    using     (6.5),    with    N   =   4      and   k  =   0,    1,     ....    63    in   the    third 

and   final    stage. 

Observe  that  in  the  second  and  third  stages  above,  the  value  of  N 
in  (6.5)  is  always  an  integer  power  of  four,  and  the  value  of  the  ex- 
ponent is  same  as  the  iteration  number.  In  total,  we  have  three 
iterations  -  computation  of  4-point,    16_point   and   64-point   DHTs  -    in  the 

P 
evaluation  of   a   64-point   DHT.      For   a   data    sequence    of    length  4    ,    we  will 

have   P   stages    in   the    evaluation  of    the   DHT  as    shown    in  Fig.    (6.1),    where 

P  =  4. 


6.4   Computational  Cost 

From  (6.5)  we  generate  the  following  seven  equations. 

H(k  +  N/4)   =  H.  (k)  -  sin(2jtk/N)n,  _  (k)   +  cos(2nk/N)H .  .,  (-k) 
4n  4n+l  4n+l 

-  cos[2n(2k)/N]Hj(   „(k)  -  sin[2n(2k)/N]H .  ^.(-k) 

4n+2  4n+2 

+   sin[2n(3k)/N]H,      _(k)   -   cos[2jt(3k)/N]H .      _(-k)(6.6) 
4n+3  4n+3 

H(k  +  N/2)      =  H,    (k)   -   cos(2jTk/N)H.      .  (k)      -    sin(2nk/N)H  (-k) 

4n  4n+l  4n+l 

+   cos[2n(2k)/N]Hj,      „(k)    +   sin[2n(2k)/N]II.    i0(-k) 
4n+2  4n+2 

-  cos[2n(3k)/N]Hv(      ,(k)   -    sin[2n(3k)/N]H        _(-k)(6.7) 

4n+3  4n+3 
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H(k  +  3N/4)  =  H.  (k)  +  sin(2nk/N)H.    (k)    -  cos(2jrk/N)H,   ,  (-k) 
4n  4n+l  4n+l 

-  cos[2n(2k)/N]H.    (k)  -  sin[2n(2k)/N]n    f-k) 

4n+2  4n+2 

-  sin[2n(3k/N]H.      _(k)    +   cos[2n(3k)/N]n        „(-k)    (6.8) 

4n+3  4n+3 

H[(N/4)   -  kj      =  H .    (-k)    +  cos(2nk/N)H.      _ (k)    +   sin(2jrk/N)H„      , (-k) 
4n  4n+l  4n+l 

+   sin[2n(2k)/N]H.         (k)    -   cos [2n(2k)/N]n,    ^(-k) 
4n+2  4n+2 

-  cos[27T(3k)/N]H.         <k)    -    sin[2n(3k)/N]ir         <-k)  (6.9) 

4n+3  4n+3 

H[(N/2)  -  k]   =  H   (-k)  +  sin(2nk/N)H.  ., (k)  -  cos (2nk/N)H,  A,  (-k) 
4n  4n+l  4n+l 

-  sin[2n(2k)/N]H.   _<k)  +  cos [2n(2k)/N]H,  ^,(-k) 

4n+2  4n+2 

+   sin[2n(3k)/N]n.      _(k)    -   cos [2*(3k)/N]H  f-k)  (6.10) 

4n+3  4n+3 

H[(3N/4)    -   k]    =   H,    (-k)    -   cos(2nk/N)Hyl         (k)    -    sin(2nk/N)R\         f-k) 
4n  4n+l  4n+l 

+    sin[2n(2k)/N]n.      .(k)    -   cos[27t(2k)/N]II .   ^.(-k) 
4n+2  4n+2 

+   cos[27T(3k)/N]H,    ^,(k)    +   sin[2jr(3k)/N]H,      ,(-k)  (6.11) 

4n+3  4n+3 

H(N  -  k)      =  H .    (-k)    +  cos(2nk/N)H.    _  (-k)   -    sin<2nk/N)H.      . (k) 
4n  4n+l  4n+l 

-   sin[2n(2k)/N]H.      _(k>    +   cos [2«(2k)/N]H .    ^-(-k) 
4n+2  4n+2 

-    sin[2n(3k)/N]H.      _<k)    +   cos[2n(3k)/N]H.    _(-k), 
4n+3  4n+3 

k  =  0,    1 (N/8)   -   1.         (6.12) 

In  deriving  the  above  equations,  we  repeatedly  used  the  fact  that 

3,  <*),  H,  ^-.(k),  Ea    ^,<k)  and  H,  ^,(k)  are  all  periodic  on  (N/4).   In 
4n      4n+l       4n+2         4n+3 

the  above  Eqns.  (6.5)  -  (6.12) 

H,  (-k)  =  H.  [(N/4)  -  k], 
4n        4n 
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V+l(-k)  ■  n4n+lt(N/2)  '  kl" 

■W-*'  "  H4n+2[<3W4)  -kl" 

H4n+3(-k)=n4n+3<N-k)- 
In  Eqns.  (6.5)  -  (6.12)  the  only  products  involved  are  those  shown  in 
Table  6.1 

TABLE  6.1  PRODUCTS  REQUIRED  TO  COMPUTE  AN  N-POINT  DUT 
USING  THE  RADIX-4  DIT  ALGORITHM 


1.   cos(27tk/N)H.   .  (k) 
4n+l 


3.   cos[2n(2k)/N]H.   _(k) 
4n+Z 

5.   cos[2n(3k)/N]H,   .(k) 
4n+3 


7.   sin(2nk/N)H,  .  (k) 
4n+l 


9.      sin[2n(2k)/N]H/(      _(k> 
4n+2 

11.      sin[2n(3k)/N]U.      _(k) 

4n+3 


2.   sin(2nk/N)H.  _ (-k) 
4n+l 


4.   sin[2n(2k)/N]n.  ^(-k) 
4n+2 

6.      sin[2n(3k)/N]IT      ,(-k) 
4n+3 

8.      cos(2nk/N)H.    _ <-k) 
4n+l 

10.   cos[2n(2k)/N]H  ^(-k) 

4n+2 

12.   cos[2n(3k)/N]H   ,(-k) 

4n+3 


k  =  0,  1,  ....  (N/8)  -  1. 
From  (6.5)  -  (6.12),  we  see  that  once  we  form  the  above  products,  we  can 
compute  the  entire  DIIT  sequence.  The  advantage  is  that  instead  of  using 
(6.5)  N  times  to  evaluate  the  N  DHT  coefficients,  we  could  as  well  use 
it  only  (N/8)  times  to  compute  the  first  (N/8)  coefficients,  and  still 
obtain  the  entire  N-point  sequence  with  the  help  of  Eqns.  (6.6)  - 
(6.12).  Table  6.2  explains  the  above  fact  for  a  64-point  sequence. 
From  Table    6.2   we    observe    that   whenever  we    compute 

H(k),      k  =   0,    1,    ...,    8      (values    in   the    left   most    column) 


using  (6.5),  with  the  help  of  the  twelve  products  in  Table  6.1,  we  can 
simultaneously  obtain  all  the  other  coefficients  in  that  row.  For 
example,  when  we  compute  R(3)  we  also  obtain  H(13),  H(19),  H(29),  H(35), 
11(45),  H(51)  and  H(61)  at  the  same  time  without  ever  actually  evaluating 
them  from  (6.5),  the  basic  DIT  equation.  This  extraordinary  symmetry  is 
not  exploited  if  we  employ  any  radix-2  algorithm  for  a  data  sequence 
whose    length   is   an   integer   power  of   four. 

TABLE      6.2      WORKING  OF  THE  RADIX-4   DIT  ALGORITHM  FOR  A  64-POINT  DHT 
k        (N/4)-k  k+N/4        (N/2)-k        k+N/2  (3N/4)-k        k+3N/4  N-k 


0 

16 

32 

48 

1 

15 

17 

31 

33 

47 

49 

63 

2 

14 

18 

30 

34 

46 

50 

62 

3 

13 

19 

29 

35 

45 

51 

61 

4 

12 

20 

28 

36 

44 

52 

60 

5 

11 

21 

27 

37 

43 

53 

59 

6 

10 

22 

26 

38 

42 

54 

58 

7 

9 

23 

25 

39 

41 

55 

57 

8 

24 

40 

56 

Therefore  we  need  12[(N/8)  +  1]  =  (3N/2)  +  2  multiplications  to 
compute  the  N-point  DHT  from  four  (N/4)-point  DHTs .  However,  whenever  k 
=  0  or  k  =  (N/8)  we  have  some  savings  in  the  number  of  multiplications. 
When  k  =  0,  there  are  no  multiplications  at  all  since  all  the 
trigonometric  functions  in  the  twelve  products  in  Table  6.1  are  either 
zero  or   one.      When  k  =    (N/8)    we   note   that 
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k  =    (N/4)    -  k  k  +    (N/4)    =    (N/2)    -   k 

k   +    (N/2)    =    (3N/4)    -   k  k   +    (3N/4)    =    (N  -   k) 

Therefore    for   this    case,    Eqns.    (6.9)    -    (6.12)    are    sane    as 

Eqns.     (6.5)    -    (6.8)     and,     as    such,    we    obtain    only    four    coefficients 

instead   of    the   usual    eight.      The    following    are    the    equations    we    obtain 

when  k     =  N/8. 

H(N/8)  =  H,  (N/8)  +   2D\   _ (N/8)  +  H.   .(N/8)  (6.13) 

4n  4n+l         4n+2 

H(3N/8)  =  H.  (N/8)  -  H.  ^„(N/8)    +   211.   .(N/8)  (6.14) 

4n         4n+2  4n+3 

H(5N/8)  =  H.  (N/8)  -   2H,   . (N/8)  +  IT  ^,(N/8)  (6.15) 

4n  4n+l         4n+2 

H(7N/8)  =  H,  (N/8)  -  H.  ^.(N/8)    -  2n.  ^.(N/8).  (6.16) 

4n         4n+2  4n+3 

In  arriving  at  the  above  equations,  we  used  the  fact 

cos(n/4)  +  sin(n/4)  =  2 

We  see  that  the  only  products  involved  in  the  above  equations  are 

2^Ea    ^-.(N/8)]  and    2[H.   .(N/8)] 
4n+l  4n+3 

So,  when  k  =  N/8,  we  need  to  compute  only  two  products  instead  of  the 
usual  twelve.  Hence  for  the  special  cases  of  k  =  0  and  k  =  N/8  com- 
bined, we  have  a  total  of  only  two  products  instead  of  the  normal  twenty 
four.   Hence  the  number  of  multiplications  M  is  given  by 

M  =  3N/2  +  12  -  22  =  (3N/2)  -  10. 
Combining  such  multiplications  over  all  the  iterations,  we  get 

M  =  [(3N/2)  -  10]  +  4[(3/2)(N/4)  -  10]  +  42 [ (3/2) (N/42)  -  10] 

+  ...  +  4P"2[(3/2)(N/4P"2)  -  10], 
where  P  =  log  N. 
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Simplifying  the  above  expression  further,  we  get 

M  =  (3N/2HP  -  1)  -  10[1  +  4  +  42  +  ...  +  4P'2] 

-  (3N/2)P  -  (3N/2)  -  10[4P_1-l]/3 
M  =  (3N/2)P  -  (7N/3)  +  10/3. 
Additions 

Once  we  form  the  twelve  products  in  Table  6.1,  let  us  form  the 
summations  given  in  Table  6.3. 

TABLE  6,3      PRIMARY  SUMMATIONS  REQUIRED  TO  COMPUTE  AN  N- POINT  DHT 

USING  THE  RADIX-4  DIT  ALGORITHM 
SUM1  =  PRODI  +  PROD2  SUM2  =  PROD3   +  PROD4 

SUM3  =  PROD5  +  PROD6  SUM4  =  PROD7   +  PROD8 

SUM5  =  PROD9  +  PRODI 0  SUM6  =  PRODI 1  +  PROD12 

where  PRODI,  PROD2,  etc.  are  the  products  in  Table  6.1.   Thus  we  have  a 

total  of  six  summations  to  compute.   Now,  again  referring  to 

Eqns.  (6.5)  -  (6.12),  we  note  that  they  can  be  implemented  as  shown 

below. 

H(k)  =  [H\  (k)  +  SUM2]   +  [SUM1  +  SUM3] 
4n 

H(k  +  N/4)    =   [H„  (k)  -  SUM2]   +  [SUM6  -  SUM4] 

4n 

H(k  +  N/2)    =   [E,  (k)  +  SUM2]   -  [SUM1  +  SUM3] 

4n 

H(k  +  3N/4)    «   [H„  (k)  -  SUM2]   -  [SUM6  -  SUM4] 

4n 

H[(N/4)  -  k]   =   [H,  (-k)  +  SUM5]  +  [SUM1  -  SUM3] 

4n 


!I[(N/2)  -  k]   =   [II.  (-k)  -  SUM5]  +  [SUM4  +  SUN6] 

4n 

R[(3N/4)  -  k]  =   [H,  (-k)  +  SUM5]  -  [SUM1  -  SUM3] 

4n 

E(N  -  k)      =   [IT  (-k)  -  SUM5]  -  [SUM4  +  SUM6] 

4n 

If  we  do  the  above  summations  directly,  we  need  three  additions  for  each 

coefficient,  thus  totaling  24  additions.   If  we  observe  the  above 

summations  carefully,  we  notice  that  each  parenthesized  sum  occurs 

twice.   So  instead  of  repeating  an  addition  already  performed,  we  store 

it  in  a  temporory  register  once  it  is  computed,  as  shown  in  Table  6.4 

TABLE  6.4   INTERMEDIATE  SUMMATIONS 

Tl  =  H,  (k)  +  SUM2  T3  =  SUM  +  SUM3 

4n 

T2  =  B.  (k)  -  SUM2  T4  =  SUM6  -  SUM4 

4n 

Then,  (6.5)  -  (6.8)  can  be  implemented  as  given  in  Table  6.5. 

TAELE  6.5  ADDITIONS  TO  GENERATE  THE  FIRST  FOUR  DKT  COEFFICIENTS 

H(k)  =  Tl  +  T3  H(k  +  N/2)  =  Tl  -  T3 

H(k  +  N/4)  =  T2  +  T4  H(k  +  3N/4)  =  T2  -  T4 

Once  H(k),  H(k  +  N/4),  H(k  +  N/2)  and  n(k  +  3N/4)  are  formed  the  sums 

stored  in  Tl  through  T4  are  no  longer  useful.   To  save  memory  space,  we 

can  use  the  same  registers  this  time  to  form  the  summations  given  in 

Table  6.6. 

TABLE  6.6      ADDITIONAL  INTERMEDIATE  SUMMATIONS 

Tl  -  E,  (-k)  +  SUM5  T3  =  SUM1  -  SUM3 

4n 

T2  =  H,  (-k)  -  SUMS  T4  =  SUM4  +  SUM6 

4n 

Then,  Eqns.  (6.9)  -  (6.12)  are  implemented  as  given  in  Table  6.7. 
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TABLE  6.7  ADDITIONS  TO  COMPUTE  THE  LAST  FOUR  DHT  COEFFICIENTS 

U[(N/4)  -  k]  =  Tl  +  T3  H[(N/2)  -  k]  =  T2  +  T4 

H[(3N/4)  -  k]  =  Tl  -  T3  H(N  -  k)  =  T2  -  T4 

If  we  add  all  the  summations  in  Tables  6.3  through  6.7,  we  have  a  total 

of  22  additions.   Normally,  we  will  have  to  carry  out  such  additions  for 

k  =  0  through  k  =  N/8   resulting  in  a  total  of 

[(N/8)  +  1]22  =  (11/4)N  +  22  additions.   However,  for  the  special  cases 

of  k  =  0  and  k  =  (N/8)  we  have  some  savings  in  the  number  of  additions. 

For  k=0,  Eqns.  (6.9)  through  (6.12)  are  only  repetitions  of  Eqns.  (6.5) 

through  (6.8)  and,  as  such,  we  obtain  only  four  distinct  coefficients. 

To  compute  them  we  form  the  four  sums  given  below. 

Tl  =  H.  (0)  +  H.  ^,(0)         T3  =  H.   .(0)  +  H.   .(0) 
4n       4n+l  4n+2       4n+3 

T2  =  H.  (0)  -  H.   .(0)         T4  =  H.    (0)  -  B.   -(0) 
4n       4n+l  4n+2       4n+3 

Then  the  four  coefficients  are  obtained  as  given  below. 

H(0)  =  Tl  +  T3  B(N/2)  =  T2  +  T4 

H(N/4)  =  Tl  -  T3  H(3N/4)  =  T2  -  4 

Thus  only  eight  additions  are  needed  when  k  =  0 

When  k  =  (N/8),  to  obtain  the  four  DHT  coefficients  given  by 

Eqns.  (6.13)  through  (6.16),  we  first  form  the  following  four 

summations. 

Tl  =  H,  (N/8)  T2  =  2[B.  ^,(N/8)] 

4n  4n+l 

T3  =  H.   0(N/8)  T4  =  2[H.  ^.(N/8)] 

4n+2  4n+3 

Then  those  equations  can  be  implemented  as  given  below 

H(N/8)  =  Tl  +  T2  +  T3 
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H(3N/8)  =  Tl  -  T3  +  T4 

H(5N/8)  =  Tl  -  T2  +  T3 

n(7N/8)  =  Tl  -  T3  -  T4. 
As  described  above,  we  thus  need  only  16  additions  for  the  special  cases 
k  =  0  and  k  =  (N/8)  combined,  instead  of  the  usual  44.   Hence  the  total 
number  of  additions  required  to  compute  the  N-point  PHT  from  four 
(N/4)-point  DKTs  is  given  by 
A  =  (11N/4)  +  22  -  44  +  16  =  (11N/4)  -  6 

Summing  all  such  additions  over  all  the  iterations  until  we  reach 
16-point  DHTs,  we  get 

A  =  [(11N/4)  -  6]  +  4[(ll/4)(N/4)  -  6]  +  42 [ (11/4) (N/42)  -  6] 

+  ...  +  4P"2[(11/4)(N/4P"2)  -  6] 
where  P  =  log  N. 

-=  ;  .   4)NP  -  1KN/4)  -  (N/2)  +  2. 

P-l 

In  addition,  we  also  have  8[4    J  =  2N  additions  to  compute  all  the 

necessary  4-point  DHTs.   Hence  the  total  number  of  additions  is 
A  =  (11N/4)P  -  5N/4  +  2. 
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7.   SPLIT-RADIX  FAST  HARTLEY  TRANSFORM  ALGORITHM 

7.1  Introduction 

In  the  previous  chapters  two  different  types  of  radix-2  algorithms 
were  described  to  compute  the  DHT  of  data  sequences  whose  length  is  an 
integer  power  of  two.  In  the  last  chapter  a  radix-4  algorithm  was 
introduced  which  computes  the  DHT  of  a  data  sequence  more  efficiently  if 
its  length  is  an  integer  power  of  four.  In  this  chapter,  a  new  algo- 
rithm which  involves  the  application  of  both  the  radix-2  and  radix-4 
algorithms  is  described.  The  new  algorithm,  known  as  the  split-radix 
algorithm,  is  computationally  more  efficient  than  either  the  radix-2  or 
the  radix-4  algorithm.  The  idea  of  the  split-radix  algorithm  was  first 
proposed  in  1984  by  Duhamel  [12]  to  compute  the  DFT  of  a  data  sequence. 
A  similar  algorithm  was  proposed  by  Soo-chang  Pei  and  Ja-ling  Wu  [13]  to 
compute  the  DHT  of  a  real-valued  sequence.  Using  the  index  mapping 
approach  Sorensen  e_t  al.  [9]  derived  a  more  efficient  split-radix  algo- 
rithm to  compute  the  DHT.  In  this  chapter  the  decomposition  formula  for 
the  algorithm  is  derived  from  the  definition  of  the  DHT  without  taking 
recourse    to   the    index  mapping   approach,    and    its   working    is   described. 

7.2  The   Decomposition  Formula 

The    split-radix   algorithm   applies    the   radix-2  DIF  decomposition   to 
the   even-indexed    samples,    and    the    radix-4    DIF   decomposition    to    the 
odd-indexed   samples   of   the   data    sequence. 
The   DHT  of   an  N-point   real-valued  data   sequence    is   given  by 


N-l 
H(k)  =  (1/N)    5  x(n)cas(27tnk/N),   k  =  0,  1,  ...,  N-l.    (7.1) 
n=0 

We  can  rewrite  the  above  summation  as  shown  below. 

(N/2)-l  N-l 

H(k)  =  (1/N)   3  x(n)cas(2nnk/N)  +  (1/N)   J  * (n)cas (2nnk/N) 
n=»0  n=N/2 

(N/2)-l  (N/2)-l 

=(1/N)[   ^  x(n)cas(2nnk/N)  +   )  x(n+N/2)cas{2nk(n  +  N/2)/N}], 
n=0  n=0 

k  =  0,  1 N-l.  (7.2) 

In  the  rest  of  the  derivation  the  factor  (1/N)  is  not  carried  through, 

but  is  restored  in  the  final  step  of  the  derivation. 

Let  us  consider  the  even-  and  odd-  indexed  DUT  coefficients 

separately. 

Let  k  =  2m  in  Eqn.  (7.2).   Then 

(N/2)-l  (N/2)-l 

H(2m)  =   ]>  x(n)cas[2:rn(2m)/N]  +   )  x(n  +  N/2)cas[2n(n  +  N/2)(2m)/N] 
n=0  n=0 

(N/2)-l  (N/2)-l 

)  x(n)cas[2nn(2m)/N]  +   ^   x(n  +  N/2)cas [2nn(2m) /N] 
n=0  n=0 

(N/2)-l 

5   fx(n)  +  x(n  +  N/2)]cas{2nnm/(N/2)}, 
n=0 

m  =  0.  1.  ....  (N/2)  -  1.   (7.3) 

The  above  is  an  (N/2)-point  DHT.   Thus  the  even-indexed  coefficients 

II(2m)  are  obtained  by  repeatedly  applying  the  above  radix-2  DIF  formula. 
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Now    let    us    consider   the   odd-indexed    samples    in   the   DllT   sequence.      Let   k 

=  2p   +   1. 

Then  we   can  rewrite  Eqn.    (7.2)    as 

(N/2)-l 
H(2p  +   1)   =       5     x(n)cas[2nn(2p  +   D/N] 
n=0 

(N/2)-l 
+        )     x(n   +   N/2)cas[2n(n  +  N/2)(2p   +   1)/N] 
n=0 

Since 

cas[(2p  +  l)n  +  x]  =  cas[n  +  x]  =  -casx, 

the  above  equation  becomes 

(N/2)-l 
H(2p  +  1)  =   ^  x(n)cas[2n(2p  +  l)n/N] 
n=0 

(N/2)-l 

^  x(n  +  N/2)cas[2n(2p  +  l)n/N]. 
n=0 

(N/2)-l 

^      [x(n)   -  x(n  +  N/2)]cas{2nn(2p  +  D/N}, 
n=0 

p  =  0,    1,    ....    (N/2)   -   1.      (7.4) 

Ilowever, 

cas[27Tn(2p  +  1)/N]   =  cos(2nn/N)cas (4nnp/N)    +   sin(2nn/N)cas{4nn(-p)} . 

Using    the    above    result    in  Eqn.    (7.4),    we    get 

(N/2)-l 
H(2p  +  1)   =       ^   [x(n)    -  x(n  +  N/2) ]cos (2nn/N)cas(4nnp/N) 
n=0 


•o  /■ 


(N/2)-l 
+   J     fx(n)  "  x(n  +  N/2)]sin(2nn/N)cas{4nn(-p)/N} . 
n=0 

If  the  direction  of  summation  in  the  second  sum  is  reversed,  we  get 

(N/2)-l 
n(2p  +  1)  =   ^   [U(n)  -  x(n  +  N/2)}cos(2jrn/N) 
n=0 

+  {x(N/2  -  n)  -  x(N  -  n)}sin(2nn/N)]cas(4nnp/N). 

P  =   0.    1,    ....(N/2)    -   1.  (7.5) 

The  above  is  the  radix-2  DIF  decomposition  formula  to  obtain  the 
odd-indexed  coefficients  of  the  DHT.  If  we  were  interested  in  only  a 
radix-2  algorithm,  we  would  have  stopped  at  this  point.  However,  in  the 
split-radix  algorithm  we  once  again  apply  a  radix-2  decomposition  to  the 
formula  in  Eqn.  (7.5).  We  observe  that  this  additional  radix-2  decom- 
position is  applied  only  in  evaluating  the  odd-indexed  coefficients,  but 
not  in  evaluating  the  even-indexed  samples.  If  we  had  done  so,  we  would 
be    effectively   deriving    the   usual    radix-4   DIF   algorithm. 

We    can    split    the     (N/2)-point    summation    in    Eqn.     (7.5)     into    two 
(N/4)-point    summations   as    follows. 

(N/4)-l 
(H(2p  +   1)   =       J  Hx(n)    -  x(n  +  N/2) }cos (2nn/N) 

n=0 

+  {x(N/2  -  n)  -  x(N  -  n)}sin(2nn/N)]cas(4nnp/N) 

(N/4)-l 
+   ^   [{x(n  +  N/4)  -  x(n  +  3N/4) }cos{2n(n  +  N/4)/N} 
n=0 

+  {x(N/4  -  n)  -  x(3N/4  -  n)}sin{2n(n  +  N/4) /N} ]cas {4np (n  +  N/4)/N} 
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Using    the    familiar    trigonometric    identities,    we    can    rewrite    the    above 

equation   as 

(N/4)-l 
n(2p  +  1)  =    J    Kx(n)  -  x(n  +  N/2)}cos(2nn/N) 
n=0 

+  {x(N/2  -  n)  -  x(N  -  n)}sin(2nn/N)]cas(4nnp/N) 

(N/4)-l 
+   5    [(x(N/4  -  n)  -  x(3N/4  -  n) }cos (2nn/N) 
n=0 

-  {x(n  +  N/4)  -  x(n  +  3N/4) } sin(2nn/N) ]cas{np  +  4nnp/N)    (7.6) 

If  we  let  p  =  2m  in  the  above  relation,  we  obtain  the  'even'  samples 

H(4m  +  1)  of  the  odd-indexed  DHT  sequence.   After  further  simplification 

we  finally  obtain 

(N/4)-l 
H(4m  +  1)  =   J    [(x(n)  -  x(n  +  N/2)  +  x(N/4  -  n)-x(3N/4-n)}cos(2nn/N) 
n=0 

+    {x(N/2-n)   -  x(N-n)    +  x(n+3N/4)    -   x(n+N/4) }  sin(27Tn/N)  ]cas  Urmia/ (N/4) } , 

m  =   0,    1,     ....    (N/4)    -   1.         (7.7) 
The   above   equation    is    the    (N/4)-point    DHT    that    computes    the    odd    DHT 
coefficients  H(4m+1),    m  =   0,    1,    ...,    (N/4)-l. 

If  we    let   p  =  2m  +   1    in   Eqn.     (7.6),    we    obtain    the    'odd'    coeffi- 
cients H(4m+3)    of   the  DHT. 

(N/4)-l 
H(4m  +   3)    =       ^  ftx(n)   -  x(n  +  N/2)}cos(2nn/N) 

nO 

+    {x(N/2   -   n)    -   x(N  -   n)}sin(2nn/N)]cas{2irn(2m  +   l)/(N/2)} 


-89- 


(N/4)-l 
+    5     [{x(n  +  3N/4)  "  l(n  +  N/4)}sin(2nn/N) 
n=0 

-  {x(N/4  -  n)  -  x(3N/4  -  n) }cos (2nn/N) ]cas {2nn(2n  +  l)/(N/2)} 

Collecting  the  terns,  we  get 

(N/4)-l 
H(4n  +  3)  =   J     {x(n)  "  *<n+N/2)  "  »(N/4  -  n)  +  x(3N/4-n) } cos (2nn/N) 
n=0 

cas{2;in(2m  +  l)/(N/2)} 

(N/4)-l 
+   ^    {x(N/2  -  n)  -  x(N  -  n)  -  x(n  +  3N/4)  +  x(n  +  N/4) } s in(2nn/N) ] 
n=0 

cas(2nn(2m  +  )/(N/2).  (7.8) 

Eowever, 

cas[2nn(2m  +  l)/(N/2)]  =  cos{2n(2n)/N}cas{2n[2nk/ (N/2) ] } 

+  sin{2n(2n)/N)cas{2n(2n)(-k)/(N/2)}.        (7.9) 

Substituting   the    above    expansion    in    the    first    sum    of    Eqn.     (7.8)    and 

calling    it    SUM1,    we    get 

(N/4)-l 
SUM1   =        2  <x(n)    "   x(n  +   N/2)    "   x(N/4   "  a)   +   x(3N/4   -   n)}cos(2nn/N) 

n=0 

cos{2n(2n)/N)cas{2n(2nk)/(N/2)} 

(N/4)-l 
+  ^  {x(n)    -   x(n   +  N/2)    -   x(N/4   -   n)    +   x(3N/4   -   n) }cos (2nn/N) 

n=0 

sin{2n(2n)/N)cas{2n(2n)(-k)/(N/2)} 
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If  we  reverse  the  direction  of  summation  of  the  second  sum  in  the  above 

equation,  i.e.  substituting  {(N/4)  -  n}  in  place  of  n,  we  get 

(N/4)-l 
SUM1  =   5    **<n>  *  s(n  +  N/2)  "  x(N/4  "  n)  +  *(3N/4  "  n)}cos(2nn/N) 
n=0 

cos[27t(2n)/N]cas{2n(2nk)/(N/2)} 

(N/4)-l 
+  5  {x(N/4   -   n)    -   x(3N/4  -   n)    -   x(n)    +  x(n  +   N/2) } s in(2nn/N) 

n=0 

sin[2n(2n)/N]  cas{2n(2nm) / (N/2) } . 

combining  the  above  two  summations  using  the  trigonometric  relation, 

cos(A  +  B)  =  cosAcosB  -  sinAsinE, 

we  get 

(N/4)-l 
SWil   =     2  <x(n)    "   x(n  +  N/2)    +  »<3N/4   -   n)    -   x(N/4  -   n) }cos{2n(3n)/N} 

n=0 

cas{2nnm/(N/4)}       (7.10) 

Now,  considering  the  second  sum  in  Eqn.  (7.8)  with  the  result  of 

Eqn.  (7.9)  substituted  in  it  and  calling  it  SUM2 ,  we  get 

(N/4)-l 
SU?!2  =   ^     {x(N/2  -  n)  -  x(N  -  n)  -  x(n  +  3N/4)  +  x(n  +  N/4)} 
n=0 

sin(2nn/N)cos{2n(2n)/N)  cas{2n(2nk)/(N/2)} 

(N/4)-l 
+   ^     {x(N/2  -  n)  -  x(N  -  n)  -  x(n  +  3N/4)  +  x(n  +  N/4)} 
n=0 
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sin(27Tn/N)sin{2n(2n)/N)    cas{2n(2n)  (-m)/(N/2)} 

reversing    the    direction    of    summation    in    the    second    sum    in    the    above 

equation,    we    get 

(N/4)-l 
SUM2   =        J  (x(N/2   -   n)    -   x(N  -   n)    -   x(n   +   3N/4)    +  x(n  +   K/4)] 

n=0 

sin(2nn/N)cos{2n(2n)/N}    cas{2n(2nx.i)  /  (N/2) } 

(N/4)-l 
+        ^  (x(N/4   +   n)    -   x(n   +   3N/4)    -   x(N  -   n)    +   x(N/2   -   n)} 

n=0 

cos(2nn/N)sin{2n(2n)/N}  cas{2n(2nn) / (N/2) } . 

Using  the  trigonometric  relation 

sin(A  +  B)  =  sinAcosB  +  cosAsinB 

and  combining  the  two  summations  in  the  above  equation  we  get 

(N/4)-l 
SUM2  =   ^      {x(N/2  -  n)  -  x(N  -  n)  -  x(n  +  3N/4)  +  x(n  +  N/4)} 
n=0 

sin(2n(3n)/N)cos{2n(nm)/(N/4)} .     (7.11) 

Combining  (7.10)  and  (7.11),  we  get 

(N/4)-l 
H(4m  +  3)  =    ^      [{x(n)  -  x(n  +  N/2)  +  x(3N/4  -  n)  -  x(N/4  -  n)} 
n=0 

cos{2n(3n)/N) 

+  {x(N/2  -  n)  -  x(N  -  n)  -  x(n  +  N/4)  +  x(n  +  N/4)} 

sin{2n(3n)/N}]cas(2nnm/(N/4), 

m  =  0,  1,  ....  (N/4)  -  1.     (7.12) 

Table  7.1  gives  a  summary  of  the  above  discussion. 
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7.3   Working  of  the  Algorithm 

Fig.  (7.1)  describes  the  working  of  the  split-radix  algorithm  in  a 
block  diagram  for  a  32-point  data  sequence.   Observe  the  way  the  DHT  is 
TABLE  7.1   PROCEDURE  TO  COMPUTE  AN  N-POINT  DHT  USING  THE  SPLIT-RADIX 

ALGORITHM 
Coefficients  to  be  computed  DIF  equation  to  be  used 

1.  Even  coefficients  (7.3) 
H(2m),  m  =  0,  1.  ....  (N/2)  -  . 

2.  Odd  coefficints  (7.7) 
H(4m  +1),   m  =  0,  1,  ....  (N/4)  -  1 

3.  Odd  coefficients  (7.12) 
H(4m  +3),   m  =  0,  1,  ...,  (N/4)  -  1. 

computed.   The  task  is  reduced  to  that  of  computing  a  half-length 
16-point  DHT  and  two  quarter-length  8-point  DHTs. 

The  half-length  DHT  computes  the  even-indexed  coefficients  H(0), 
H(2).  H(4),  H(6),  H(8),  H(10),  H(12),  H(14),  H(16),  H(18).  H(20).  H(22). 
H(24),  H(26),  n(28),  H(30) . 

The  first  quarter-length  DHT  computes  the  odd-indexed  coefficients 
H(l).  E(5),  H(9),  H(13),  H(17).  H(21),  H(25),  H(29) . 

The  second  quarter-length  DHT  computes  the  following  odd-indexed 
coefficients  H(3).  H(7),  n(ll),  H(15),  H(19),  H(23),  H(27),  H(31). 
However  the  computation  of  the  even-  and  odd-  indexed  coefficients  is 
again  reduced  to  three  smaller  length  DHT  computations  until  we  reach 
the  stage  where  we  need  to  compute  only  4-point  or  2-point  DHTs,  both  of 


-93- 


which  require  no  multiplications.  The  split-radix  algorithm  also  is  an 
in-place  algorithm.  Tt  is  applicable  to  all  sequences  whose  length  is 
an  integer  power  of  two,  and  is  the  most  efficient  of  all  in-place 
algorithms. 
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The  number  of  omit ipl icat ions  N  is  given  by 

M  =  (2N/3)log2N  -  19(N/9)  +  3  +  (-l)P/9 
and  the  the  number  of  additions  A  is  given  by 

A  =  (4N/3)log2N  -  (14N/9)  +  3  +  (-l)P(5/9) 

where  P  =  log.N. 

Tables  7.2  through  7.4  give  a  summary  of  the  operation  counts  for 
different  algorithms.  By  operation  counts  we  mean  the  total  number  of 
multiplications  and  additions.  We  can  see  from  these  operation  counts 
that  the  split-radix  algorithm  uses  less  multiplications  and  less  addi- 
tions than  either  the  radix-2  or  the  radix-4,  but  still  has  a  fairly 
compact  code.  Therefore  this  is  generally  the  most  efficient  algorithm 
for  data  sequences  whose  length  is  a  power  of  two.  The  operation  counts 
for  these  algorithms  is  less  than  those  required  to  compute  the  DFT  of  a 
real-valued  data  sequence  [14], 
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TABLE  7.2  OPERATION-COUNTS  FOR  RADIX-2  ALGORITHMS 


Length 

Mults 

Adds 

Mult  +  Add 

RADIX-2 

4 

0 

8 

8 

8 

4 

26 

30 

16 

20 

74 

94 

32 

68 

194 

262 

64 

196 

482 

678 

128 

516 

1154 

1670 

256 

1284 

2690 

3974 

512 

3076 

6146 

9222 

1024 

7172 

13826 

20998 

2048 

16388 

30722 

47110 

4096 

36868 

67586 

104454 

TABLE  7.3   OPERATION-COUNTS  FOR  THE  RADIX-4  ALGORITHM 


Length 

Mults 

Adds 

Mult  +  Add 

4 

0 

8 

8 

16 

14 

70 

84 

64 

142 

450 

594 

256 

942 

2498 

3440 

1024 

5294 

12802 

18096 

4096 

27310 

62466 

89776 
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TABLE  7.4  OPERATION-COUNTS  FOR  THE  SPLIT-RADIX  ALGORITHM 


Lenp.th 

Mults 

Adds 

Mult  +  Add 

SPLIT-RADIX 

4 

0 

8 

8 

8 

2 

22 

24 

16 

12 

64 

76 

32 

42 

166 

208 

64 

124 

416 

540 

128 

330 

998 

1328 

256 

828 

2336 

3164 

512 

1994 

5350 

7344 

1024 

4668 

12064 

16732 

2048 

10698 

26854 

37552 

4096 

24124 

59168 

83292 
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8.   AN  APPLICATION:   RAMAN  SPECTRA  AND  MATCHED  FILTERING 

8.1  Introduction 

In  this  chapter  an  example  application  of  the  DHT  in  the  analysis 
of  Raman  spectra  is  described,  wherein  it  is  employed  to  perform  the 
fast  convolution  operation.  We  begin  the  chapter  with  a  brief  introduc- 
tion to  the  origin  of  Raman  spectra,  and  then  we  discuss  the  role  of  the 
FT  and  the  DFT  in  the  resolution  of  the  spectra.  The  concept  of  matched 
filtering  and  its  applicability  in  the  enhancement  of  SNR  of  the  spectra 
is  presented.  The  implementation  of  the  matched  filter  using  both  the 
DFT  and  the  DHT  is  discussed  next,  and  finally  we  close  the  chapter  with 
a  comparison  of  the  performance  of  the  DFT  and  the  DHT  in  the  applica- 
tion chosen. 

8.2  Raman  Spectra 

The  Raman  effect,  also  known  as  the  scattering  of  modified  radia- 
tion, was  first  discovered  in  1928  by  C.  V.  Raman  [15].  The  phenomenon, 
first  observed  in  the  case  of  liquids,  is  universal  in  character  and  is 
a  powerful  technique  in  identifying  different  substances  unambiguously. 
When  visible  light  of  a  particular  frequency  is  used  to  excite  any 
material,  the  spectrum  of  the  scattered  light  from  the  substance  con- 
tains the  signature  of  the  material.  This  phenomenon  is  known  as  the 
Raman  effect,  and  the  spectrum  of  the  scattered  light  is  referred  to  as 
the  Raman  spectrum  for  that  substance.  For  any  given  substance,  the 
difference  in  frequency  between  the  incident  light  and  the  scattered 
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light  is  constant  and  is  independent  of  the  frequency  of  the  incident 
light.  This  shift  in  frequency  between  the  incident  light  and  the 
scattered  light  corresponds  to  the  frequency  of  oscillation  of  the 
chemically  bonded  atoms  in  the  substance,  which  in  turn  depends  on  the 
geometry  of  the  molecule.  In  terms  of  quantum  theory  of  light,  the 
difference  in  frequency  between  the  spectrum  of  the  incident  light  and 
the  Raman  spectrum  is  positive  or  negative,  depending  on  whether  the 
incident  light  is  delivering  energy  to  the  target  molecule  or  receiving 
energy  from  it. 

One  major  drawback  in  analyzing  the  Raman  spectra  is  that  they  are 
very  feeble.  Matched  filtering  is  a  very  useful  technique  to  enhance 
the  SNR  in  such  cases. 

S.3   Matched  Filtering 

8.3.1  Preliminaries 

Before  describing  the  matched  filter,  some  preliminary  definitions 
and  models  for  the  signal  and  noise  portion  of  the  spectrometer  output 
are   presented. 

Raman  shift:  When  a  substance  is  subjected  to  the  Raman  effect,  the 
frequency  of  oscillation  of  the  chemically  bonded  atoms  in  the  substance 
alters.   This  shift  in  frequency  is  called  the  Raman  shift  and  is 

measured  in  terms  of  wavenumbers  (cm   ) . 
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signals:  Spectroscopic  signals  are  pulse-like  in  nature  with  wavenumber 
as  the  independent  variable.  For  our  study  here  a  Lorentzian  pulse  is 
considered.      It    is   defined  by 

s(v)    =     A(vQ)/    [1   +    {(v  -  vQ)/h}2] 

where  A(v  )  is  the  peak  amplitude  of  s(v)  at  the  center  frequency  v  and 

h  is  the  half-width  at  half-height  (hwhh)  of  the  pulse.   The  Lorentzian 

pulse  is  taken  for  study  since  it  is  the  most  difficult  of  all  the  peaks 

to  detect  for  a  given  half-width  at  half-height. 

Noise:   The  noise  part  of  the  spectrum  is  treated  as  a  random  process 

(rp)  with  white  Gaussian  statistics. 

White  noise:   So  named  by  analogy  to  the  white  light  which  contains  all 

the  visible  frequencies,  white  noise  is  a  zero  mean  rp  whose  power 

spectral  density  (PSD)  is  a  constant  at  all  the  frequencies. 

PSD:   If  x(t)  is  a  wide-sense  stationary  (wss)  rp,  and  its  FT  is  given 

by  X(f),  then  its  PSD  S   (f)  is  given  by 


xz 


XX 


T  — >  » 


E  |  |XT(f)!2  /  2T  ]. 


where  7^(f)  is  the  FT  of  the  signal  x(t)  restricted  to  an  interval  -T 

and  T.   For  the  white  noise  S   (f)  is  a  constant  and  is  given  by 

xx 

S   (f)  =  NA  /  2.  (8.1) 

xx       0 

Autocorrelation:       If    x(t)     is    a   wss    rp,     its    autocorrelation   function 
R      (t)    is    given  by 
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R 
xx 


(t)   -  FT*1   [   S      (f)    ].  (8.2) 


From  Eqns.    (8.1)    and    (8.2)    we    see    that    R       (t)     for    the    white    noise     is 

xx 


given  by 


R   (t)  =  [N  /2]6(t), 
XX  o 


where,  6(t)  is  the  Dirac  delta  function.   As  is  clear  from  the  above 
equation,  the  white  noise  is  uncorrelated  to  itself. 

8.3.2  Theory  of  the  Matched  Filter 

Matched  filter  is  an  optimal  linear  system  which  maximizes  the  SNR 
at  the  output  of  the  filter. 

We  assume  that  the  input  to  the  filter  consists  of  the  sum  of  a 

deterministic  signal  s(t)  and  a  random  noise  process  n(t).   The  signal 

and  noise  at  the  output  of  the  filter  are  denoted  as  s  (t)  and  n  (t), 

o  o 

respectively.      The    instantaneous    signal    power   at    the    output    is    given  by 

2 
s      (t)    and    the    average    noise    power    as    N    .       If    the    filter's    impulse 
o  o 

response    is    indicated    as    h(t)    and    its   Fourier    transform   as    11(f),    then 

the    output    signal    s    (T)    at    t   =  T    is    given  by 


S(f)H(f)ej2,TfT  df. 


-  I 


where  S(f)  is  the  FT  of  the  signal  s(t). 
Similarly  the  output  average  noise  power  is 
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Nq   =      J     S^mllim^df  (8.3) 

—CO 

where    S.,.,(f)     is    the    PSD    of    the    input    noise    process.      Therefore    the 

NN 

output    SNR   is 

00 

S   /N     =      |     f   S(f)H(f)ej27TfT  df    I2    /   N 
o     o  J  o 


where   N      is    given    in  Eqn.    (8.3) 
o 

It    can  be    shown    [16]    that    the    above   ratio    is  maximum  when 

•m  -J2nfT 

E(f)    =  H        (£)    =   C   &   ).',-. (8.4) 

opt  S^(f) 

where  C  is  a  real  constant  and  S  (f)  is  the  complex  conjugate  of  S(f). 


8.3.3   Matched  Filter  for  White  Noise 

We  know  that  if  the  input  noise  process  is  white,  its  PSD 

S»,»,(f)  =  N  12.      Substituting  the  above  value  for  SXIXt(f)  in  Eqn.  (8.4), 

NN       o  NN 

we  get 

H  At)   =  KS*(f)e"j2TTfT 
opt 

where  K  =  2C/N   is  a  real  constant.   The  impulse  response  of  the  filter 
o 

is  given  by 


h    ft)  =  f  H    (f)ej2nft  df 
opt      J    opt 
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=  Ks(T  -  t).  (£.5) 

From  Eqn.  (8.5)  it  is  clear  that  the  impulse  response  of  the  matched 

filter,  when  the  input  noise  process  is  white,  is  equal  to  the  input 

signal  shifted  by  T  units  and  then  reversed.   The  delay  T  is  a  parameter 

under  the  control  of  the  filter  designer.   If  we  choose  T  =  0  and  K  =  1 

in  Eqn.  (8.5),  then  it  becomes 

h    ft)  =  s(-t).  (8.6) 

opt 

Since  the  transfer  function  of  the  filter  is  so  closely  tied  to  the 

input  signal  as  shown  in  the  above  equation,  the  filter  is  referred  to 

as  the  matched  filter.   This  feature  can  sometimes  be  a  problem  too,  in 

that  it  requires  of  us  a  knowledge  of  the  exact  shape  of  the  input 

signal.   However,  we  can  often  make  good  guesses  about  the  input  signal 

profile . 

8.3.4   Matched  Filter  as  a  Correlator 

The  output  of  the  matched  filter  is  given  by 

y(t)  =  x(t)*h    ft) 
opt 


=  f  x(r)h    ft  -  x)dx 
J      opt 


where  x(t)  =  s(t)  +  n(t)  and  *  denotes  convolution.   If  n(t)  is  a  white 
noise  process,  the  above  equation  becomes 


y(t)  =  J 


x(x)s(-t  +  r)  dt 
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»  J  x(t   +  T)s(r)dT  (8.7) 

—00 

The  above  equation  shows  that  y(t)  is  the  cross  correlation  of  the 
assumed    input    signal  with   the    input    signal-plus-noise. 

8.3.5   Implementation  Considerations 

The  spectrum  to  be  filtered  is  treated  as  an  N-point  sequence 

{x(m)},     a   =  0,  1,  ...,  N  -  1. 

and  the  sampled  version  {h    (m)}  of  the  filter's  impulse  response  is 

opt 

also  made  the  same  length.   Then  the  linear  convolution  of  the  input  and 

the  filter's  impulse  response  is  given  by 

N-l 
y(k)  =  3  x<m)h(k  -  m).  k  =  0,  1,  ....  2N  -  2,  (8.8) 

m=0 

where  y(k)  is  the  output  of  the  matched  filter.  The  convolution  equa- 
tion above  is  the  linear  convolution  of  the  input  spectrum  and  the 
filter's  impulse  response.  However,  the  linear  convolution  can  be 
achieved  by  adding  atleast  N-l  zeros  to  both  {x(m)}  and  {h(m)}  and 
then  performing  a  circular  convolution  on  the  resulting  sequences.  The 
advantage  of  the  approach  lies  in  the  fact  that  we  could  use  the  DFT  or 
the  DHT  to  compute  the  circular  convolution.  Fig.  (8.1)  shows  the  block 
diagram  for  implementing  the  matched  filter  using  the  DFT.  Dyer  and 
Hardin  [17]  used  the  above  approach  in  analyzing  simulated  Raman  spectra 
using  the  matched  filter  technique. 
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Implementation  of  the  matched  filter  using  the  DI'T  approach  is 
presented   here. 

Fig.  (8.2)  shows  the  block  diagram  of  the  implementation  of  fast 
linear  convolution  using  the  DKT.  The  N-point  spectrum  {x(m)}  and  the 
N-point  digitized  filter  transfer  function  {h(m)}  are  both  zero  padded 
upto  2N  points.  Mote  that  for  the  matched  filter,  {h(m)}  is  the 
reversed  version  of  the  input  pulse.  Since  a  Lorentzian  pulse  centered 
at  zero  is  an  even  pulse,  obeying  the  relation  h(m)  =  h(N  -  m),  the  Din 
of  the  circular  convolution  of  x(m)  and  h(m)  is  the  product  of  the 
individual    DKTs    H    (k)    and   H    (k)    of   x(m)    and   h(m),    respectively.      Using 

this  fact,  we  find  the  DIIT  of  the  circular  convolution  of  the  2N-point 
sequences  (x(m)}  and  (h(m)}  as  shown  in  Fig.  3.2.  Since  the  DIIT  is  a 
symmetric  transform,  by  performing  the  2N-point  DIIT  on  the  resulting 
sequence,  we  get  the  linear  convolution  of  x(m)  and  h(m) .  Since  we  are 
interested  only  in  the  enhancement  of  the  N-point  spectrum  that  we 
started  with,  we  discard  the  last  N  points  in  the  resulting  sequence. 
The  main  advantages  of  the  DIIT  approach  are  that  it  computes  the  linear 
convolution  faster  than  the  DFT  technique  and,  that  it  is  simpler  to 
implement  than  the  latter  when  one  of  the  functions  being  convolved  is 
even. 
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8.4      Simulation  Results 

As    an   example    application   of   the   DIIT,    simulated  Raman   spectra  were 
analyzed   using   matched   filter.    The    signal    taken   for    study    consisted    of 

two   Lorentzian   peaks,    one    centered    at    250    cm        and   the   other  at   650 

cm      ,    and  was    shown    in  Fig.    8.3.      The   half-width   at    half-height    (hwhh) 

of  each  peak  was  20  cm  ,  and  the  sampling  rate  v/as  1024  sanples/sec . 
Yfhite  Gaussian  noise  having  unity  variance  was  added  to  the  signal  to 
generate  the  input  spectrum.  Fig.  8.4  represents  the  noise  and  Fig.  8.5 
depicts   the    input    spectrum. 

The  input  spectrum  consisted  of  1024  points  and  the  SNR  (maximum 
amplitude  of  peak/  rms  amplitude  of  noise)  was  2.0.  Fig.  $.6  presents 
the  output  of  a  matched  filter  having  the  spectrum  of  Fig.  8.5.  as  its 
input.      The    impulse   response    of    the    filter   was    a    reversed   Lorentzian 

peak    having    a    hwhh    of    20    cm       .      The    two    spectral   peaks   are    evident    in 

the   output    spectrum,    one    located   at   250   cm       and   the   other   at    650    cm 

Fig.  8.7  presents  the  output  of  the  same  matched  filter  iznlemented 
using    the   DIIT   instead   of   the   DFT.      Again,    the    spectral   peaks    are    clearly 

discernible    from   the   noise    at   250   cm        and   650   cm 

As  evident  from  the  output  spectrum,  the  filter  output  is  insensi- 
tive to  the  choice  of  the  transform.  The  analysis  was  carried  out  with 
two   other   filters   also,    one    matched    to    a    Lorentzian    peak    of    hwhh   =    5 
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cm  and  the  other  to  a  Lorentzian  peak  with  hwhh  =  10  cm  .  The  fil- 
ters v/ere  implemented  with  the  DFT  as  well  as  the  D!1T.  Table  8.1  shows 
a  summary  of  the  computation  times  for  both  the  methods. 

TABLE  8.1   COMPUTATION  TIME  FOR  IMPLEMENTING  THE  MATCHED  FILTER 

HWHH  of  filter's  response    Approach       CPU  time  in  sees. 

<-   -1 
5  cm 


10  cm"1 


20  cm 


D!IT 

2.80 

DFT 

3.06 

DIIT 

2.83 

DFT 

3.00 

DHT 

2.77 

DFT 

3.05 

On  the   average    the    DHT   approach   has    taken   2.8    sees    compared    to    3.04 
sees,    for   the   DFT   approach   for   a   1024   point    spectrum. 

Software  for  the  matched  filter  was  written  in  VAX-11  FORTRAN  and 
was  given  in  the  Appendix.  The  white  noise  was  generated  using  RALPH,  a 
general-purpose    DSP    software    developed    at    Kansas    State    University. 
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9.   CONCLUSIONS 

Various  properties  of  the  HT  for  a  real  signal,  and  the  DIIT  for  a 
real  data  sequence  were  derived  in  the  first  two  chapters.  It  was  shown 
that  all  the  properties  of  the  FT  and  the  DFT  have  a  counterpart  in  the 
HT  theory.  These  properties  look  very  much  identical,  especially  when 
the  data  sequence  possesses  either  odd  or  even  symmetry. 

The  DFT  of  a  data  sequence  generates  a  complex-valued  sequence. 
Ilowever,  the  DFT  of  a  real-valued  data  has  a  useful  property  that  its 
real  part  is  even  and  its  imaginary  part  is  odd.  General-purpose  FFT 
algorithms  written  to  deal  with  real-  as  well  as  complex-valued  data  do 
not  take  advantage  of  the  above  property.  However  those  algorithms  can 
be  optimized  for  computing  the  DFT  of  real-valued  data  sequences  more 
efficiently  to  save  memory  and  computation  time.  The  DHT  helps  us  avoid 
taking  recourse  to  such  complicated  methods  of  optimizing  to  achieve 
efficiency  in  dealing  with  real-valued  data.  Once  the  DHT  of  the 
real-valued  data  is  computed,  using  a  simple  relation  between  the  DFT 
and  the  DHT,  we  obtain  the  DFT  of  the  sequence.  This  approach  for 
computing  the  DFT  has  the  advantage  of  speedier  computation  of  the  DIIT, 
and  it  does  not  have  the  complexity  of  optimizing  methods.  This  method 
is  especially  attractive  in  case  of  a  DIF  FFT  for  real-valued  input, 
since  there  are  no  easy  methods  available  to  compute  the  DFT  of  a  real 
data  efficiently  from  the  conventional  methods.  The  operation-counts 
for  radix-4  and  split-radix  algorithms  to  compute  the  DHT  also  are  less 
than  the  number  of  operations  for  the  corresponding  FFT  algorithms. 
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Finally  an  esample  application  was  given,  in  which  simulated  Raman 
spectra  were  analyzed  using  a  matched  filter.  The  DKT  approach  took 
less      computational      effort      than     the     DFT     method. 
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APPENDIX 


ft****************  t-******^*******************************,,,****** 

Department  of  Electrical  and  Computer  Engineering 

Kansas  State  University 
Vax  Fortran  Source  Filename:  MATCHED_FILTER.FOR 

*******e**«*«*******************************«*****************«« 
ROUTINE:      Mainline 

DESCRIPTION:  It  enhances  the  Signal  to  Noise  Ratio 
(SNR)  of  the  input  spectrum  using  the 
Matched  Filtering  method. 
Discrete  Hartley  Transform  (DHT) 
technique  is  used  to  implement  the  matched 
filter 

DOCUMENTATION 
FILES :        None 

RETURN:       Not  Used 

ROUTINES 

CALLED:  DHT2,NPLOT 

AUTHOR:        CHANDRA  C.  VARANASI 

DATE  CREATED:  22nd  March  1987,   Version  1.0 

***«**«*************«************>M:  *************************** 

The  following  are  some  of  the  key  variables  used  in 
the  routine 

SPECTRUM:  (Real)  Array  containing  the  input  spectrum, 
zero  padded  to  twice  its  original 
length 

PROFILE:   (Real)  Array  containing  the  signal  profile 

pulse, zero  padded  to  twice  its  original 
length 

OUTPUT:    (Real)  Array  containing  the  spectrum  with 

enhanced  Signal  to  Noise  Ratio  (SNR), 
resulting  as  the  output  from  the 
matched  filter. Length  of  the  array  is 
same  as  the  original  input  spectral 
data  length 
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IMPLICIT  NONE 

REM.  SPECTRUM (0 : 2047 ). PROFILE (0: 2047 ) , X_DATA (0:2047), 
+       OUTPUT(0:1023),Y_DATA(0:2047) 

INTEGER  NUMBER, I 

*  Read  the  input  spectral  data  into  the  array  SPECTRUM 

CALL  SGOPEN  ( 5  ,  '  READ  ' ,  '  NO  PROMPT ' ,  '  SPECT .  DAT ' ,  '  REAL ' .  NUMBER ) 
CALL  SGTRAN  ( 5 ,  ' READ ' , ' REAL ' . SPECTRUM. NUMBER ) 

DO  I  =  0,2047 

Y_DATA ( I )  =  SPECTRUM ( I ) 

X_DATA(I)=FLOAT(I) 
END  DO 

*  Read  the  input  signal  profile  pulse  into  the  array  PROFILE 

CALL  SGOPEN  (6, 'READ' , 'NOPROMPT' , 'PRO. DAT' , 'REAL' , NUMBER) 
CALL  SGTRAN  ( 6 ,  ' READ ' , ' REAL ' , PROF ILE , NUMBER ) 

*  Take  the  DHT  of  the  SPECTRUM 

CALL  LIB  INIT_TIMER 

CALL  DHT2  (SPECTRUM, 2048, 11 ) 

*  Take  the  DHT  of  the  signal  PROFILE 

CALL  DHT2  (PROFILE, 2048. 11) 

*  Multiply  the  above  two  DHTs 

DO   I  =  0,2047 

SPECTRUM(I)    =    SPECTRUM (I )*PROFILE( I) 
END  DO 

*  Take    the    inverse    DHT  of    the    sequence. In   fact, inverse   DHT  and 

*  forward  DKT  are   one   and   the    same 

CALL  DHT2    (SPECTRUM, 2048,11 ) 
CALL  LIB   SnOW_TI?IER 

*  Now, retain  the   first    1024   points 

DO   I  =  0,1023 

OUTPUT(I)    =    SPECTRUM(I) 

X_DATA(I)    =  FLOAT (I) 
END  DO 


-122- 


*  Plot  the  matched  filter  OUTPUT 

CALL  S IMPLK_PLOT (1024.' 
+  'LINEAR ',X_DATA. OUTPUT. 'cm  I t-lt I '. 

+  'Raman  Shift', '  '.'Amplitude') 

END 
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*«*********(l*******#«***«*******£4******4:«  *********************** 

Department  of  Electrical  and  Computer  Engineering 

Kansas  State  University 
Vax  Fortran  Source  Filename:  MATCDED_FILTER_DFT.FOR 

t************  ********** +  +:***.>:■*(■**  +  ****•**  +  +  *  +  •*  +  <!'****  *********** 
ROUTINE:       Mainline 

DESCRIPTION:   It  enhances  the  Signal  to  Noise  Ratio 
(SNR)  of  the  input  spectrum  using  the 
Matched  Filtering  method. 
Discrete  Fourier  Transform  (DFT)  is  used 
to  implement  the  matched  filter 

DOCUJENTATION 
FILES:        None 

RETURN:       Not  Used 

ROUTINES 

CALLED:        FFT.NPLOT 

AUTHOR:        CHANDRA  C.  VARANASI 

DATE  CREATED:  22nd  March  1987.     Version  1.0 

************************************************************** 
The  following  are  some  of  the  key  variables  used  in 
the  routine 

SPECTRUM:  (Complex)  Array  containing  the  input  spectrum, 

(zero  padded  to  twice  its  original 
length),  as  its  real  part . Imaginary 
part  is  zero. 

PROFILE:   (Complex)  Array  containing  the  signal  profile 

pulse, (zero  padded  to  twice  its 
original  length),  as  its  real  part. 
Imaginary  part  is  zero 

OUTPUT:    (Real)     Array  containing  the  spectrum  v/ith 

enhanced  Signal  to  Noise  Ratio  (SNR), 
resulting  as  the  output  from  the 
matched  filter. Length  of  the  array  is 
same  as  the  original  spectral  length 

*****************************  *><****f.***>:*K  ***********  ********** 
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IMPLICIT  NONE 

COMPLEX   SPECTRUM* 0 : 2047 ) .PROFILE (0 : 2047) 

REAL  OUTPUT(0:1023),X_DATA(0:1023),A(0:2047),B(0:2047) 

INTEGER   NUMBER, I 

CALL  LIB  INIT_TIMER 

*  Read  the  input  spectral  data  into  the  real  array  A 

CALL  SGOPEN ( 5 ,  ' READ ' » ' NOPROMPT ' , ' SPECT . DAT ' , ' REAL ' , NUMBER ) 
CALL  SGTRAN  (5. 'READ' , 'REAL' , A, NUMBER) 

*  Read  the  zero  padded  profile  pulse  into  the  real  array  B 

CALL  SGOPEN  ( 6 .  ' READ ' , ' NOPROMPT ' , ' PRO . DAT ' , ' REAL ' , NUMBER ) 
CALL  SGTRAN  (6. 'READ' , 'REAL' ,B, NUMBER) 

*  Transfer  the  data  and  profile  into  the  complex  arrays  SPECTRUM 

*  and  PROFILE  respectively 

DO  I  =  0,2047 

PROFILE ( I )=  B(I) 

SPECTRUM(I)  =  A(I) 
END  DO 

*  Take  the  DFT  of  the  input  SPECTRUM 

CALL  LIB  INIT_TIMER 

CALL  FFT  (SPECTRUM, 2048,0) 

*  Take  the  conjugate  of  the  transformed  SPECTRUM 

DO  I  =  0,2047 

SPECTRUM(I)  =  CONJG(SPECTRUM(I)) 
END  DO 

*  Take  the  DFT  of  the  signal  PROFILE 

CALL  FFT  (PROFILE, 204 8,0) 

*  Multiply  X  *(K)  and  H(K)  where: 

*  X(K)   DFT  of  the  zero  padded  inpuf  spectrum 

*  IKK)   DFT  of  the  zero  padded  signal  profile  pulse 

DO  I  =  0,2047 

SPECTRUM(I)  =  SPECTRUM (I) 'PROFILE (I) 
END  DO 
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*  Take  the  inverse  DFT  of  the  product 

CALL  FFT  (SPECTRUM, 2048,1) 
CALL  LIB  SE0',7_TIMER 

*  Now  retain  only  1024  points  of  the  data,  since  that  is  the 

*  length  of  the  original  spectrum 

DO  I  =  1,1023 

OUTPUT(I)  =REAL(SPECTRUM(204S-I)) 

X_DATA(I)  =  REAL(I) 
END  DO 

*  Plot  the  matched  filter  output 

CALL  SIHPLE_PLOT( 1024, 'Matched  Filter  Output  using  DFT', 
+  ' LINEAR ',X_D ATA, OUTPUT. 'cm|t-lt| ', 

+  'Raman  Shift','  ' , 'Ampl itude ' ) 

END 
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*  4c  *  ft  ft  ft  *  *  *  *  *  *  ft  ft  ft  ft  s*  *******  *  *  *  *  ft  ft  *  *  *  *  *  *  ft  ft  ft  ft  ft  ft  *  *  ft  *  *  *  *  ij  *************  * 

Department  of  Electrical  and  Computer  Engineering 

Kansas  State  University 
Vax  Fortran  Source  Filename:   PROFILE. FOR 

************************** ft*  ft*  *  *  *  *  >s  ft  *  *  *  *  *  >s  *  *  »s  *  ft  ft  >s  >:••  ft  e  ft  *  c  ft  ******** 
ROUTINE:      PROFILE. FOR 

DESCRIPTION:   It  generates  a  Lorentzian  pulse  centered 
at  zero.   After  that,  the  array 
containing  the  pulse  is  zero  padded 
to  twice  the  pulse  length.   Zeros  are 
added  in  the  middle  of  the  data  strean 
instead  of  at  the  end.   The  zero  padded 
pulse  will  be  used  as  the  profile  in  the 
convolution 

DOCUMENTATION 
FILES:        None 


RETURN: 

ROUTINES 
CALLED: 


Not  Used 


None 

AUTHOR:       CHANDRA  C.  VARANASI 
DATE  CREATED:  30th  March  1987    Version  1.0 

**************  *  * ***********  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  ft  ft  *  *  *  *  *********  * 

The  following  are  the  key  variables  used  in  the  routine 
AMPLITUDE: 


ALPHA: 

PULSE: 
PROFILE : 


(Real)   The  amplitude  of  the  Lorentzian 
peak 

(Real)   The  half_width  at  half_height 
of  the  pulse 

(Real)   Array  containing  the  pulse 

(Real)   Array  containing  the  zero  padded 
pulse 


******************  *  *  ft  *  ft  ft  ft  ft  ft  **********  *  *  *  ft  ft  ft  ft  ft ft  ft  ft  ft  ft  ft  ft  C ft  ft  ft  *  *******  « 
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IMPLICIT  NONE 

REAL  AMPLITUDE,  ALPIIA,  PULSE(-512 : 511 ) ,  PROFILE(0:  2047  )  , 
+    X_DATA( 0:2047) 

INTEGER    I 

PARAMETER  (AMPLITUDE  =  2.0.  ALPHA  =  20.0) 

*  Generate  the  Lorentzian  pulse 

DO  I  =  -512,511 

PULSE(I)  =  AMPLITUDE/ (1+  (I/ALPHA)**2) 
END  DO 

*  Now,  zero  pad  the  profile  upto  twice  the  number  of  points 

*  bringing  the  negative  part  of  the  pulse  to  the  end  of  the 

*  data  stream  and  store  in  the  array  PROFILE 

DO  I  =  -512,-1 

PROFILE (2048+1)  =  PULSE(I) 

PROFILE (1+5 12)   =  PULSE (1+5 12) 
END  DO 

*  Write  PROFILE  to  the  disk 

CALL  SGOPEN  (5, 'WRITE ' , 'NOPROMPT' , 'PRO. DAT ' , 'REAL' , 2048) 
CALL  SGTRAN  ( 5 , ' WRITE ' , ' REAL ' »  PROFILE , 204  8 ) 

*  Plot  the  zero_padded  profile.   For  that  set  up  the  X-axis  data 

DO  I  =  0,2047 

X_DATA(I)  =  REAL(I) 
END  DO 

*  Now  plot 

CALL  SIMPLE_PLOT(2048, 'Zero  Padded  Prof ile ', 'LINEAR' , 
+  X_DATA, PROFILE, 'cm|t-lt I ', 'Raman  Shift', 

+  '  ', 'Amplitude') 

END 
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************  *  *  *  «  *  #  *  *  *  *  *  *  **£►:.  *  ****.;:►:;  * *  *  *  *  *  *  *  *  a  *  *  *  *  #  *  *  *  #  *  *  *  *  >:,  $ ^.  *  ** 

Department  of  Electrical  and  Computer  Engineering 

Kansas  State  University 
Vax  Fortran  Source  Filename:   NOISE_PLOT.FOR 

ROUTINE :      NOI SE_PLOT . FOR 

DESCRIPTION:   It  plots  the  zero  mean  unit  variance 
white  Gaussian  noise  generated  using 
RALPH,  a  general  purpose  signal 
processing  software  developed  at  the 
Kansas  State  University. 

DOCUMENTATION 
FILES :        None 

RETURN:       Not  Used 

ROUTINES 

CALLED :       S IMPLE_PLOT 

AUTHOR:       CHANDRA  C.  VARANASI 

DATE  CREATED:  30th  March  1987    Version  1.0 

*********************************************************** * *« ** 

The  following  are  the  key  variables  used  in  the  routine 

GAUSSIAN  NOISE:     (Real)Array  containing  the  samples  of 

the  zero  mean  unit  variance  Gaussian 
noise . 

*  ♦  *  *  *  *  *  *  *  *  *  *  *  *  lie  A  *  *  *  ft  ft  ft  *  *  *  ft  ft  *  *  *  *  ft  *  *ft  *  ft  ft*  *  ft*  *  *  *  *  *  ft*  *  ft  ft  *  *  ft  ft  ft  ft  *  *  *  *  *  ft 

IMPLICIT  NONE 

REAL       X_DATA (0:1023).  GAUSSIAN  NOISE (0: 1023 ) 

INTEGER    I,  NUMBER 

PARAMETER   (AMPLITUDE  =2.0,  ALPHA  =  20.0) 

*  Read  the  noise  samples  from  the  disk 

CALL  SGOPEN  ( 5 ,  ' READ ' , ' NOPROMPT ' , ' NOI SE . DAT ' , ' REAL ' , 
+  NUMBER) 
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CALL  SGTRAN  ( 5 ,  ' READ ' , ' REAL ' . G AU SS I AN_NO I SE . NUMBER ) 

DO  I  =  0,  1023 

X_DATA(I)  =  FLOAT(I) 
END  DO 

*  Plot  the  noise 

CALL  SIMPLE_PLOT( 1024, 'White  Gaussian  noise ', 'LINEAR' , 
+  X_DATA,GAUSSIAN_NOISE,  'Cm|t-lt|'» 

+  'Raman  shift',  '  ',  'Amplitude') 

END 
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Department  of  Electrical  and  Computer  Engineering 

Kansas  State  University 
Vax  Fortran  Source  Filename:  SPECTRUM. FOR 

**************  **4^***«***«********«*«**««*****«*«$$« «&«&$#««**«« 

ROUTINE:      SPECTRUM. FOR 

DESCRIPTION:   It  reads  the  input  signal  data  file  as 
well  as  the  noise  data  file  and  then 
generates  the  input  spectrum  by  nixing 
the  above  two  files 

DOCUMENTATION 
FILES :        None 

RETURN:       Not  Used 

ROUTINES 

CALLED :  S IMPLE_PLOT 

AUTHOR:       CHANDRA  C.  VARANASI 

DATE  CREATED:  1st  April  1987     Version  1.0 

***********«****«****#*#**********************«******«********** 

The  following  are  the  key  variables  used  in  the  routine 

SPECTRUM:      (Real)   Array  into  which  the  input  signal 

is  read.   When  the  noise  is  added 
to  it,  it  contains  the  simulated 
input  spectrum  to  the  matched 
filter 

NOISE:        (Real)   Array  into  which  the  white 

Gaussian  noise,  generated  using 
RALPH,  is  read 

*<««*iti4(*a******#*<:*<:H:<:«:**««*«************************************** 

IMPLICIT    NONE 

REAL  SPECTRUM(0:2047),NOISE(0:2047),X_DATA(0:2047) 
INTEGER    NUMBER,  I 
*  Read  the  input  signal  into  SPECTRUM.   Input  signal  is  a 
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*  misture  of  two  Lorentzian  pulses. 

CALL  SGOPEN  ( 7 ,  ' READ ' , ' NOPROMPT ' , ' S IGNAL . DAT ' , ' REAL ' , 

+  NUMBER) 

CALL  SGTRAM  ( 7 ,  ' READ ' , ' REAL ' , SPECTRUM , NUMBER ) 

*  Read  the  zero  mean  unit  variance  white  Gaussian  noise  generated 

*  using  RALPH 

CALL  SGOPEN  (5, 'READ' , 'NOPROMPT' , 'NOISE.DAT '. 'REAL' , 

+  NUMBER) 

CALL  SGTRAN  ( 5 , ' READ ' , ' REAL ' . NO I SE , NUMBER ) 

*  Add  the  two  arrays  SPECTRUM  and  NOISE  to  generate  the  input 

*  spectrum  and  place  the  result  in  SPECTRUM 

DO  I  =  0.2047 

SPECTRUM (I)  =  SPECTRUM(I)+NOISE(I) 

X_DATA(I)  =  REAL(I) 
END  DO 

*  Write  the  input  spectrum  to  the  disk 

CALL  SGOPEN  (6, 'WRITE ', 'NOPROMPT' , 'SPECT.DAT' , 'REAL' ,2048) 
CALL  SGTRAN  ( 6 ,  ' WRITE ' , ' REAL ' , SPECTRUM ,2048) 

*  Plot  the  input  spectrum 

CALL  SIMPLE_PLOT( 2048, 'Input  Spectrum ', 'LINEAR' ,X_DATA. 

+        SPECTRUM, 'cm  It-It  I ', 'Raman  Shift', '  '.'Amplitude') 


END 
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*  *  *****  *  ************************** ******************************* 


Department  of  Electrical  and  Computer  Engineering 
Kansas  State  University 
Vax  Fortran  Source  Filename:      LORENTZ.FOR 

**************************************************************** 


ROUTINE: 
DESCRIPTION: 


DOCU!!ENTATION 
FILES: 

RETURN: 

ROUTINES 
CALLED: 

AUTHOR: 

DATE  CREATED: 


LORENTZ.FOR 

It  generates  a  lorentzian  shaped  pulse  of 
desired  amplitude  and  half  width  half 
height 


None 
Not  Used 

None 

CHANDRA  C.VARANASI 

23rd  March,  1987    Version  1.0 


**************************************************************** 

The  following  are  some  of  the  key  variables  used  in  the 
routine 


PEAK: 


ALPHA: 


(Real)  The  amplitude  of  the 
lorentzian  peak 

(Real)  The  half_width  at  half_height 


CENTER_FREQUENCY :  (Real)  The  frequency  at  which  the 

lorentzian  pulse  attains  its 
peak 

L0RENTZ1,L0RENTZ2: (Real)  Arrays  that  contain  the 

lorentzian  pulses  at  two 
distinct  center  frequencies 

LORENTZ:  (Real)  Array  that  contains  the 

summation  of  L0RENTZ1  and 
LORENTZ2 

************************>:  *************************************** 


1*%  * 


I'fPLICIT  NONE 


REAL 


+ 
+ 
+ 
+ 


PEAK,  ALPHA,  CENTER_FREQUENCY1 , 
CENTER_FREQUENCY2 , LORENTZ1 (0 : 1023 ) , 
FREQUENCIES ( 0 : 1023 ) ,  LORENTZ2 ( 0 : 1023 ) , 
LORENTZ (0:1023). 
X  DATA(0:1023) 


INTEGER 
PARAMETER 


+ 
+ 


(PEAK  =  2.0,  ALPHA  =  20.0, 
CENTER_FREQUENCY1  =  250.0, 
CENTER  FREQUENCY2  =  650.0) 


*  *  **********************  *  *  * »::  *  $  *  >: a .  >::  <•  c  ■:  •  >:  >:•  >:•  *  *  *  *  *  *  *  *  *  **************  * 


DO  I  =  0.1023 

LORENTZl(I)  =  PEAK/(1.0+((I  -  CENTER_FREQUENCY1 ) / 

ALPHA) **2) 
LORENTZ2(I)  =  PEAK/ (1.0+( (I  -  CENTER_FR£QUENCY2 ) / 

ALPHA) **2) 


END  DO 

*  Add  LORENTZ1  and  LORENTZ2   that  is  going  to  be  the  signal 

DO  I  =  0,1023 

LORENTZ (I)  =  L0RENTZ1(I)+L0RENTZ2(I) 

X_DATA(I)   =  REAL(I) 

FREQUENCIES (I)  =  REAL(I) 
END  DO 

*  Plot  the  signal 

CALL  SIMPLE_PLOT( 1024, 'Input  Signal ', 'LINEAR  * ,X_D ATA, 
+  LORENTZ, 'cmlt-ltl '» 'Raman  Shift', '  '» 

+  'Amplitude') 

*  Write  the  signal  to  the  disk 

CALL  SGOPEN  ( 6 , ' WRITE ' , ' NOPROMPT ' , ' S IGNAL . DAT ' . ' REAL ' , 
+  1024) 

CALL  SGTRAN  (6, 'WRITE '» 'REAL' .LORENTZ, 1024) 

END 
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#  ft  *****  ft  ft  *  ft  ft  ft  ft  ft  *  ft  ft  *  ft  ft  ft  ft  ft  *  ft  >:=  ft  *  ********►::**  ft  ft  ft  ft  ft  *******************  * 
ft 

Department  of  Electrical  and  Computer  Engineering 

Kansas  State  University 
Vax  Fortran  source  filename:      SIMPLE_PLOT.FOR 

ft*************************************************************** 

ROUTINE:  SUBROUTINE 

SIMPLE_PLOT  (NUM_POINTS.PL0T_TrTLE, 
PLOTJYPE ,  X_DATA»  Y_DATA,  X_AXI  S_UNITS , 
X_AXI S_TITLE , Y_AXI SJJNITS , Y_AXI S_TITLE ) 

DESCRIPTION:     Plots  X_DATA  and  Y_DATA  values  as 

abscissae  and  ordinates  respectively,  and 
places  a  TITLE  and  UNITS  on  the  corresp- 
onding axes.   PLOT_TYPE  specifies  the 
type  of  plot  being  generated, and  finally 
the  plot  is  given  a  PLOT_TITLE. 
The  routine  gives  the  user  the  choice  of 
plotting  on  either  the  screen  or  the  HP 
plotter. 

DOCUMENTATION 
FILES:  None 

VARIABLES  IN 
THE  ARGUMENT: 


NUM  POINTS: 


PLOT  TITLE: 


PLOT  TYPE: 


(input)   integer 

Number  of  data  points  to  be  plotted 

(input)   character*(*) 
Title  to  be  placed  on  plot 

(input)   character*(*) 

Character  string  specifying  the  type  of 

plot  to  be  generated.   The  following  are 

valid: 

for    linear-linear 

for  log-linear 

for  linear-log 

for  log-log 


'LINEAR' 
'LOG-LINEAR' 
'LINEAR-LOG' 
'LOG-LOG' 


X  DATA: 


Y  DATA: 


(input)    real 

Array  of  abscissae  values  to  be  plotted 

(input)    real 

Array  of  ordinate  values  to  be  plotted 
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X_AXIS_UNITS:  (input)    character* (* ) 

Name  to  be  given  to  the  units  associated 
with  the  X_AXIS 

X_AXIS_TITLE:  (input)    character* (* ) 

Title  to  be  placed  on  the  X_AXIS 

Y_AXIS_UNITS:  (input)    character* (*) 

Name  to  be  given  to  units  associated  with 
Y_axis 

Y_AXIS_TITLE:  (input)    character* (* ) 

Title  to  be  placed  on  the  Y_axis 


RETURN: 

ROUTINES 
CALLED : 


Not    used 


PAXIS 

PCHRPL 

PCLOSP 

PINIT 

PLGAXS 

PLGLIN 

PLGLOG 

PLINE 

PLOGSC 

PL^OG 

PORIG 

PPLOT 

PSCALE 

PSTCIIR 

PSTVEL 

PTEXT 

PTXTLN 

PWIND 


AUTHOR: 

DATE   CREATED: 

REVISIONS: 


Chandra  C.  Varanasi 
2nd  September  1986 


8th  September  1986 

Plotting  log-log  curve  is  added 

13th  September  1986 

Offsetting  the  plot_title  is  added 

**»«*«t*pii********«*#«»*««»iii#*i5*ii(*(tf)i(«^o*f:t:.^*i::i:i:ft:c^f(>:'<<i****«** 

The  subroutine  SIMPLE_PLOT  calls  some  routines  from 
THE  P  SYSTEM  OF  GENERALIZED  PLOT  SUBROUTINES.   The 
following  are  the  key  variables  required  to  call  the 
above  routines. 
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CLENX: 


CLENY: 


(Output  from  subroutine  PLOGSC) 
Length  in  cm.  of  one  log  cycle 
1.0  <=  CLEN  <=  Length  of  X_axis 

(Output  from  subroutine  PLOGSC) 
Length  in  cm.  of  one  log  cycle 
1.0  <=  CLEN  <=  Length  of  Y_Axis 


real 


DELTAX: 


(output  from  subroutine  PSCALE)    real 
Scale  factor  (increment  in  X  between  tic 
marks) 


DELTAY: 


(output  from  the  subroutine  PSCALE)  real 
Scale  factor  (increment  in  Y  between  tic 
marks) 


DIVLENX: 


(output  from  subroutine  PSCALE)     real 
Space  between  X_azis  tic  marks  in  user 
units 


DIVLENY: 


(output  from  subroutine  PSCALE)  real 
Space  between  Y_axis  tic  marks  in  user 
units 


FIRDEL: 


(input  to  suj .  -  line  PLINE)  real 
A  four  element  array  containing  the 
following: 


FIRSTX 
DELTAX 
FIRSTY 
DELTAY 


Starting  value  of  X_data 
Described  above 
Starting  value  of  Y 
Described  above 


FIRLEN: 


(input  to  subroutine  PLINE)     real 

A  four  element  containing  the  following: 


FIRSTX 
DELTAX 
FIRSTY 
CLENY: 


Starting  value  of  X_DATA 
Described  above 
Starting  value  of  Y_DATA 
Described  above 


LENSTR: 


(output  from  subroutine  PTXTLN)    integer 
Number  of  characters  in  PLOT  TITLE 


NEGFLGX : 


(output  from  subroutine  PLOGSC)   integer 

Flag  to  warn  that  negative  values  were 

encountered  in  X_DATA 

0:   No  negative  values  in  X_DATA 

1:   Negative  values  in  X_DATA 
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NEGFLGY ! 


OFFSET: 


(output  from  subroutine  PLOGSC)   integer 
Flag  to  warn  that  negative  values  were 
encountered  in  Y_DATA 

(input  to  the  routine  PPLOT)    real 
Margin  on  the  left  side  of  the  title) 


«***«*******««*****<^*************«4t*****************  ********** 


SUBROUTINE  SIMPLE_PLOT (NUM_POINTS, PLOT_TITLE, 

+  PLOTJTYPE , X_D ATA , Y_DAT A , 

+  X_AXI S_UNITS . X_AXI S_TITLE . 

+  Y_AXI S_UNTTS ,  Y_AXI  S_TITLF. ) 

Declare  the  variables 

IMPLICIT  NONE 

INTEGER  CHOICE. I , NUM_POINTS, NEGFLGX, NEGFLGY. 

+  LENSTR 


REAL  CLENX,  CLENY,  DELTAX.  DELTAY,  DIVLENX. 

DIVLENY.  FIRDEL(4).  FIRLEN(4).  FIRSTX. 
FIRSTY.  OFFSET.  X_DATA(0:*).  Y_DATA(0:*) 

CHARACTER     PLOTJITLE*  (*) ,  PLOTJTYPE* ( * ) , 

X_AXIS_UNTTS* (* ) . X_AXIS_TITLE* (• ) 
Y  AXIS  UNITS*(*),  Y  AXIS  TITLE*(*) 


**************************************************************** 

PRINT*.  '  SELECT  THE  PLOTTING  DEVICE  (ENTER  1  OR  2 ) ' 

PRINT* 

PRINT*,  '  1.   Tektronix  4014  display' 

PRINT*,'  2.   HP  7475A  Plotter' 

*  Read  the  user's  choice  of  the  plotting  device 

READ*.  CHOICE 

*  Initiate  the  corresponding  device 

IF  (CHOICE  .EG.  1)  THEN 

CALL  PINTT (4014, '  ',1.0, 'A') 
ELSE 

CALL  PINIT( 7475, '  ',1.0, 'A') 
END  IF 
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*  Set   pen  velocity 

CALL  PSTVEL    (3.0) 

*  Establish  an  origin 

CALL  PORIG  (4.54,  4.0) 

*  Establish  the  bounds  of  the  plot 

CALL  PWIND  (0.0,0.0,0.0,0.0) 

IF  (PLOTJTYPE  .EQ.  'LINEAR')  THEN 

*  Scale  both  X_DATA  and  Y_DATA 

CALL  PSCALE  (X_DATA, NUM_POINTS ,18.0, FIRSTX. DELTAX, 
+  DIVLENX) 

CALL  PSCALE  (Y_DATA, NUM_POINTS, 12 . 0 , FIRSTY , DELTAY, 
+  DIVLENY) 

*  Fill  in  the  array  FIRDEL 

FIRDEL(l)  =  FIRSTX 
FIRDEL (2)  =  BELTAX 
FIRDEL(3)  =  FIRSTY 

FIRDEL (4)  =  DELTAY 

* 

*  Draw  the  linear  axes 

CALL  PAXIS(0.0,0.0,X_AXIS_TITLE,X_AXIS_UNITS,220. 
2010, 18. 0,0.0, FIRSTX, DELTAX, DIVLENX) 

CALL  PAXIS(0.0,0.0,Y_AXIS_TITLE,Y_AXIS_UNITS,120, 
1010, 12. 0,90.0, FIRSTY, DELTAY, DIVLENY) 

*  Draw  the  linear-linear  curve 

CALL  PLINE  (X_DATA, Y_DATA» NUM_POINTS, FIRDEL, 1 , '  ' , 
+  DIVLENX, DIVLENY) 

END  IF 

IF  (PL0T_TYPE  .EQ.  'LOG-LINEAR')  THEN 

*  Scale  the  X  and  Y  axes  data 

CALL  PSCALE (X_DATA , NUM_POINTS ,18.0, FIRSTX , DELTAX , 
+  DIVLENX) 
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CALL  PLOGSC(Y_DATA,  NUMJPOINTS,  12 .0, FIRSTY,  CLEW , NEGFLGY ) 

*  Draw  the  logarithmic  Y  and  linear  X  axes 

CALL  PAXIS (0 .0,0.0, X_AXIS_TITLE, X_AXIS_UNITS, 220 , 2010 , 
+  18. 0.0.0, FIRSTX. DELTAX.DIVLENX) 

CALL  PLGAXS(0.0,0.0,Y_AXIS_TITLE,Y_AXIS_UNITS,-1010, 
+  12. 0,90.0, FIRSTY, CLENY) 

*  Fill  in  the  array  FIRLEN 

FIRDEL(l)  =  FIRSTX 

FIRDEL(2)  =  DELTAX 

FIRDELO)  =  FIRSTY 

FIRDEL(4)  =  CLENY 

*  Draw  the  log-linear  curve 

CALL  PLGLIN (X_DATA, Y_DATA , NUM_POINTS , FIRLEN, 1 ,  ' 
+  DIVLENX) 

END  IF 


IF  (PLOT_TYPE  .EQ.  'LINEAR-LOG')  THEN 

*  Scale  the  data 

CALL  PLOGSC(Y_DATA, NUM_POINTS. 18 .0, FIRSTX, CLENX, NEGFLGX) 
CALL  PSCALE(Y_DATA, ND?.l_POINTS,  12  .0, FIRSTY, DELTAY, 
+  DIVLENY,  ) 


*  Draw  the  axes 

CALL  PAXIS  (0.0,0.0,Y_AXIS_TITIJE,Y_AXIS_UNITS,  120,1010, 
+  12. 0,90.0, FIRSTY, DELTAY,D I VLENY) 

CALL  PLGAXS(0.0,0.0,X_AXIS_TITLE,X_AXIS_UNITS, 2010, 18.0, 
+  0.0, FIRSTX, CLENX) 

* 

*  Fill  in  the  array  FIRLEN 

FIRDEL(l)  =  FIRSTX 

FIRDEL(2)  =  CLENX 

FIRDELO)  =  FIRSTY 

FIRDELU)  =  DELTAY 

*  Draw  the  linear-log  curve 

CALL  PLNLOG (X_DATA, Y_DATA, NUM_POINTS , F JRLEN, 1 , ' 
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+  DIVLENY) 

END  IF 

IF  (PLOTJTYPE  .EQ.  'LOG-LINEAR')  THEN 

* 

*  Scale  the  data 

CALL  PLOGSC (X_DATA, NUM_POINTS, 18 .0, FIRSTX, CLENX, NEGFLGX) 
CALL  PLOGSC (YX_DATA.NUM_POINTS, 12 .O.FIRSTY, CLENY , NEGFLGY ) 

*  Draw  the  logarithmic  axes 

CALL  PLGAXS(0.0,0.0,X_AXIS_TITLE,X_AXIS_UNITS, 2010, 18.0, 
+  0.0, FIRSTX, CLENX) 

CALL  PLGAXS (0 . 0, 0 . 0 , Y_AXIS_TITLE, Y_AXIS_UNITS, -1010,12 . 0 , 
+  90.  O.FIRSTY,  CLEW) 

*  Fill  in  the  array  FIRLEN 

FIRDEL(l)  =  FIRSTX 

FIRDEL(2)  =  CLENX 

FIRDELO)  =  FIRSTY 

FIRDEL(4)  =  DELTAY 

*  Plot  the  log-log  curve 

CALL  PLGLOG (X_DATA, Y_DATA, NUM_PO I NTS. FIRLEN, 1 , '  ' ) 
END  IF 

*  Calculate  the  length  of  the  plot_title 

CALL  PTXTLN  (PLOT_TITLE,LENSTR) 

*  Set  the  character  size 

CALL  PSTCER  (0.3,0.4,10.0) 

*  Set  the  offset 

OFFSET  =  (18.0  -  LENSTR*0.3)1.5)/2.0 

*  Move  the  pen  to  the  position  where  the  plot_title  is  to  start 

CALL  PPLOT(OFFSET,  13.0,  0) 

*  Plot  the  title 

CALL  PTEXT(PLOT  TITLE) 
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*  Close  the  plotting  device 

CALL  PCLOSP 

END 
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ABSTRACT 

The  discrete  Fourier  transform  (DFT)  is  widely  used  in  various 
applications  in  digital  signal  procesing  (DSP).  However,  the  DFT  is  a 
complex-valued  transform  and,  as  such,  it  transforms  real-valued  se- 
quences into  complex-valued  transform  sequences.  Also,  the  inverse  DFT 
is  different  from  the  forward  DFT.  In  contrast,  the  discrete  Hartley 
transform  (DnT)  is  an  inherently  real-valued  transform  and  is  also 
symmetric . 

Various  properties  of  the  Hartley  transform  (IIT)  and  the  DHT  are 
derived  in  this  study  and  compared  with  the  well-known  properties  of  the 
Fourier  transform  (FT)  and  the  DFT.  Decomposition  formulas  for  dif- 
ferent fast  algorithms  to  compute  the  DHT  are  derived.  The  algorithms 
include  the  decimation-in-time  (DIT)  and  decimat  ion- in-f  requency  (DIF) 
radix-2,  radix-4  and  split-radix  algorithms.  Computational  cost  in 
terms  of  number  of  multiplications  and  additions  is  derived  for  all  the 
algorithms  and  compared.  Finally,  an  application  of  the  DHT  in  the 
analysis  of  Raman  spectra  is  given,  wherein  it  is  used  to  implement  the 
matched  filter.  The  matched  filter  was  implemented  via  the  DFT  also. 
Both  filter  response  and  computation  time  were  compared  for  the  two 
methods.  The  DHT  approach  proved  faster  in  terms  of  computation  time, 
but    the    filter   response   was    insensitive    to   the    choice    of   the    transform. 


